R 小筆記 – 2

R 真是不一般, 首先它不像別的語言, index 是從 0 開始. 而是從 1 開始. 再來呢, 很多地方的 1 看起來是 1 但不是 1; 看起來是 0 但也不是 0; 2 也不是平常的 2. 底下舉例說明.

上次說到用 X predict Y, 像是身高預測. 因為最低身高不是 0, 這條 minimum error 的線是有截距的. 此時一般線性迴歸的斜率公式是:

  • 模型:y = α + βx
  • 斜率:β = r(x,y) × sd(y)/sd(x)

假設強制要求通過原點, 他的迴歸如下.

  • 模型:y = βx(無截距項)
  • 斜率:β = Σ(xy)/Σ(x²)

最後, 假設原本是有截距的迴歸, 但我把它的數據先做正規化, 那它會跟誰一樣嗎?

x <- c(1.2, 2.8, 3.5, 1.9, 4.1, 2.3, 3.8, 1.6, 2.9, 3.2)
y <- c(2.4, 5.1, 6.8, 3.2, 7.9, 4.1, 7.2, 2.8, 5.5, 6.1)

# 標準化(z-score normalization)
x_norm <- (x - mean(x)) / sd(x)
y_norm <- (y - mean(y)) / sd(y)

# 方法0:一般線性迴歸斜率公式
slope0 <- cor(x, y) * sd(y) / sd(x)

# 方法1:通過原點迴歸
slope1 <- sum(x * y) / sum(x^2)

# 方法2:先正規化後當作一般線性迴歸
slope2 <- cor(x_norm, y_norm) * sd(y_norm) / sd(x_norm) 

cat("有截距 - 不通過原點迴歸:", slope0, "\n")
cat("通過原點迴歸:", slope1, "\n")
cat("標準化後做一般回歸:", slope2, "\n")

結果我們會得到 slope0 ≈ 1.85, slope1 ≈ 1.94, slope2 ≈ 0.93, 三個都不一樣.

方法 0 和方法 2 當中, 因為按照定義 :

sd(x_norm) = sd(y_norm) = 1, mean = 0

而 correlation 方面: cor(x,y) 又等於 cor(x_norm, y_norm),

所以 sd(y)/sd(x) ≠ 1 的情況下, 會導致 slope0 ≠ slope2.

方法 1 是強制通過原點, 而原本最佳解可能是有截距的, 所以 slope0 也不太容易巧合和 slope1 、slope2 一樣.

上面講的還沒有發揮到 R 的精神, 我們用 lm() 重做一次.

x <- c(1.2, 2.8, 3.5, 1.9, 4.1, 2.3, 3.8, 1.6, 2.9, 3.2)
y <- c(2.4, 5.1, 6.8, 3.2, 7.9, 4.1, 7.2, 2.8, 5.5, 6.1)

# 標準化
x_norm <- (x - mean(x)) / sd(x)
y_norm <- (y - mean(y)) / sd(y)

# === 手動計算(原始方法) ===
slope0_manual <- cor(x, y) * sd(y) / sd(x)
slope1_manual <- sum(x * y) / sum(x^2)
slope2_manual <- cor(x_norm, y_norm) * sd(y_norm) / sd(x_norm)

# === 使用 lm() 實現 ===
# 方法0:一般線性迴歸(有截距)
model0 <- lm(y ~ x)
slope0_lm <- coef(model0)[2]  # 取斜率係數

# 方法1:通過原點迴歸(無截距)
model1 <- lm(y ~ x - 1)  # -1 表示移除截距項
slope1_lm <- coef(model1)[1]  # 只有一個係數

# 方法2:標準化後的一般迴歸
model2 <- lm(y_norm ~ x_norm)
slope2_lm <- coef(model2)[2]  # 取斜率係數

# === 比較結果 ===
cat("=== 手動計算 vs lm() 比較 ===\n")
cat("方法0 - 有截距迴歸:\n")
cat("  手動計算:", slope0_manual, "\n")
cat("  lm() 結果:", slope0_lm, "\n")
cat("  差異:", abs(slope0_manual - slope0_lm), "\n\n")

cat("方法1 - 通過原點迴歸:\n")
cat("  手動計算:", slope1_manual, "\n")
cat("  lm() 結果:", slope1_lm, "\n")
cat("  差異:", abs(slope1_manual - slope1_lm), "\n\n")

cat("方法2 - 標準化後迴歸:\n")
cat("  手動計算:", slope2_manual, "\n")
cat("  lm() 結果:", slope2_lm, "\n")
cat("  差異:", abs(slope2_manual - slope2_lm), "\n\n")

# === 顯示完整模型資訊 ===
cat("=== 完整模型摘要 ===\n")
cat("\n--- 方法0:一般線性迴歸 ---\n")
print(summary(model0))

cat("\n--- 方法1:通過原點迴歸 ---\n")
print(summary(model1))

cat("\n--- 方法2:標準化後迴歸 ---\n")
print(summary(model2))

我們可以看到 model1 <- lm(y ~ x – 1) , 它不是真的減 1, -1 表示移除截距項. 還可以表示為 y ~ 0 + x. 兩者等效.

再來 slope2_lm <- coef(model2)[2] 的 [2] 不是第三個係數. 其中, 第一個係數 [1] 就固定是截距, 第二個係數 [2] 就固定是斜率.

Rrrrrrrrrrrrrrrr 啊啊啊啊啊啊啊啊啊啊….. (聲音漸尖)

R 小筆記

Realtek 公司也簡稱 R, 不過應該沒有多少人用過 R 語言. 最近適逢要用到 regression model, 但網路課程幾乎都是用 R (歷史包袱), 而不是 Python, 所以我也試著瞭解用 R 來理解 regression model.

首先 library(UsingR) 載入基本套件, 然後用 x 來 predict y.

y = β0 + β1 * X

在程式中, 以 beta0 取代 β0 , 而 beta1 取代β 1 , 因為我們只是 predict, 並不是 assign 值給 y, 所以要求的是 minimum square error.

i=1n​(Yi−(β0​+β1Xi))2

為了增加趣味性, 課程舉例是用父母的身高預測小孩的身高. 我們可以理解到身高不會為 0, 所以有一個 beta0 撐著, 然後父母身高 X 乘上的 weighting beta1, 就可以預測小孩的身高. 當然, 我們又可以理解到這是在統計上才有意義, 對個案沒有意義.

lm(y~x) 就可以得到用 x 預測 y 的 linear model 的結果.

coef(lm(y~x)) 回報 coefficients beta0 和 beta1.

c(beta0, beta1) 用來產生向量.

舉例來說, x 和 y 用 c() 產生向量, 得出 linear model 後, model 再用 coef() 取出係數, 然後印出來.


x <- c(1, 2, 3, 4, 5)
y <- c(2, 3, 5, 7, 11)

# 建立線性回歸模型
model <- lm(y ~ x)

# 獲取模型係數
coefficients <- coef(model)  # This gives c(beta0, beta1)

intercept <- coefficients[1] # beta0
slope <- coefficients[2] # beta1

# 印模型係數
print(intercept)
print(slope)

# OR 一行印模型係數
print(coefficients)

其中 <- 表示 assign, 正式名稱為賦值運算子 (assignment operator), 和 C 的 pointer 反方向.

用 Python 其實也不會變囉嗦. 以這個例子來說, Python 可以輕鬆應對. 雖然 R 確實更簡潔一點點.

import numpy as np
from sklearn.linear_model import LinearRegression

x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 3, 5, 7, 11])

# 建立線性回歸模型
model = LinearRegression()
model.fit(x.reshape(-1, 1), y)

# 獲取模型係數
intercept = model.intercept_
slope = model.coef_[0]

# 印模型係數
print(f"Intercept (beta0): {intercept}")
print(f"Slope (beta1): {slope}")

# OR 一行印模型係數
print(model)
 

美國國庫券小數學 – 2

上次寫那篇主要是因為 AI 出錯, 引發了我想記錄下來的興趣. 但後來看到一些類題需要腦筋再多轉幾個彎. 我還是做個筆記才好複習.

Q1, 有個 T-Bills 用 99,900 買了, 還有 180 天到期. yield to maturity 是多少?

上回講到計算隱含價值, 這回我買都買了, 狀況就不一樣. 所以不是用那個公式, 倒是可以用直覺.

由於 180 天後我可以拿到 100,000 元, 表示我可以賺 100,000 – 99,900 = 100 (元).

獲利率 = 100/99,900, 後續可以當作 Annual Payment Rate = APR 來計算.

年化 yield to maturity = 100/99,900 x (360/180) = 0.2002%

Q2. EAR 是多少?

Effective Annual Rate = EAR = (1 + APR/m)m , 本例子中, 我剛好買在生命周期的一半, 故兩個半年相當於兩期, m = 2.

EAR = (1 + 0.1001%)360/180 – 1 = 0.2003%

這明顯是一個買貴的例子. 但注意 EAR 只會比 yield to maturity 大一點點. 如果數量級差太多就是算錯了.

Q3. 隔夜逆回購 (ON RRP = Over Night Reverse Repurchase)

A 銀行在美聯儲的存款帳戶(reserve account)有閒置的 10 億美金, 向美聯儲申請參與 ON RRP, 假設隔夜逆回購利率 5.8%, 1 天能賺到多少錢?

美聯儲接受 ON RRP 申請後, 銀行將指示美聯儲從自己在聯儲的帳戶扣款 (減少 reserve balance), 當作國債的抵押品. 銀行則得到國債一紙.

第二天, 美聯儲把錢還給銀行, 並且支付國債利息  10億美元 * 5.80% * (1/360) ≒ 16.11 萬美元. A 銀行把國債還給美聯儲.

假設匯率 1 USD = 29.49, TWD, 16.11 萬美元 = 475.1 萬元台幣.

參考 [1], 銀行不勞而獲的結果, 剛好可以養一個螃蟹工程師 (對, 數字我湊的). 我的重點是, 錢每天都在貶值, 你怎麼可以不投資?

[REF]

  1. https://news.tvbs.com.tw/life/2889863

美國國庫券小數學

美國國庫券 (Treasury Bills, T-Bills) 價錢怎麼算呢? 其實很簡單, 但是要計算天數. 我本來想把這個惱人的事情給 AI 做, 想不到最近數學大有進步的 AI 竟然就瘋了?

[ 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 + 30 + 31 + 31 + 28 + 30 = 3 + 31 + 30 + 31 +…. 後面塞爆自動停了.

我只好自己把上面加起來, 以利息來算, 經過 244 天.

Bo = FaceVelue * (1 – d * D / 360), 其中

FaceValue = 面值

d= 折扣率 = discount quote (%),

D = 到期日 = day to mature (貨幣市場都小於一年)

360 = 約定俗成的一年天數

Bo = 隱含價格

Discount quote 是什麼呢? 根據報價網站 [1],我們可以看到 US Treasury Quotes 的數據. 記得要選到 Treasure Bills (下圖紅色畫底線), 另外一邊是 Treasury Notes.

其中, bid (買價) 和 asked (賣價) 都是 discount quote. 所以會算出兩個隱含價值. 在網路課程中, 舉的 discount 都是 0.4%, 過了七八年, discount 已經變成 4% 了!

利率變高好像沒有什麼大不了的, 債券利率也很明顯比以前變高啊! But, 債券是有債息的, 國庫券根本不配息, 它只有面值. 只會透過低於面值的折扣價來買賣.

舉例來說, 美聯儲發行 Treasure Bills 讓銀行買, 銀行用折扣價買入, 準備賺個年利率 5%. 當銀行現金不足, 可以把 T-Bills 賣給別的銀行周轉. 這個流動性是債券所沒有的, 因為它一年之內就會到期, 等同自帶 timer 的現金.

若是經濟不好的時候, 美國不用發現金、不用發 5 倍劵, 而是加價買回市場上的 T-Bills. 因為銀行就是賺錢機器, 在利益優先的考慮之下, 就自動成了美聯儲的打手, 喔不, 推手. 當銀行變得滿手現金, 自然就會想要用它賺更多錢, 也就間接刺激了市場成長. 當然這是最理想的狀況. 等我學成之後, 應該可以舉出很多歪樓的狀況.

[REF]

  1. https://www.wsj.com/market-data/bonds/treasuries

買賣債券小數學

purchased a 15-year bond that pays semi-annual coupon of $20 and is currently selling at par. What would your realized annual return be if you sold the bond five years later when the yield is 5.5%?

這題看起來很短, 但是足夠複習一些知識了.

原本大家都想知道美國降息之後, 債券會漲多少? 這題剛好相反, 是問上升後, 逢高賣出可賺多少? 我們要先知道賣出的合理價錢. 也就是 5 年後的 Pv . 會讓我腦筋打結的地方是要從 15 年後回推 5 年後, 所以是 10 年期間, 半年配息一次, 共 20 個週期. 假設面值 1,000.

Pv = 20 * ADF(5.5%/2, 20) + 1,000 / (1+5.5%/2)^20

上回筆記到 Annual Discount Factor = ADF(r, n) = (1−(1+r)-n)/r

因此在利率上漲到 5.5% 時, 該債券值 885.7956. 此處假設未來第 6~15 年利率都是 5.5%.

因為是按照面值買的 (at par), 所以我的成本是 1,000. 最後價值是 885.8, 求年報酬率 g, 半年是 g/2.

1000 = 20 * ADF(g/2, 10) + 885.7956 / (1+g/2)^10

沒有財務計算機就用 Excel. 其中 905.8 = 885.8 + 20, 最後一年是本利合計.

g/2 =IRR([-1000, 20, 20, 20, 20, 20, 20, 20, 20, 20, 905.80])

解出 g = 1.807%.


知道這個技巧後, 換個場景.

20 年美債, 面值 1000 元, 殖利率 5%, 半年配, 以現值 Pv = ADF(5%/2, 40) + 1,000 / (1+5%/2)^40 = 1000 合理價買入.

買了持有 5 年, 利率果真降到 3%, 此時債券價值上漲為 1240.158 . 公式跟上面都一樣.

Pv = 1240.158 = 20 * ADF(3%/2, 30) + 1,000 / (1+3%/2)^30

這段期間的 IRR([-1000.39, 25, 25, 25, 25, 25, 25, 25, 25, 25, 1240.158+25]) = 4.45%

整年 = 半年 * 2 = 8.9%.

假如債券殖利率真的這樣發展, 跟每年台灣大盤的平均漲幅也差不多. 但抱得愈久, 或是殖利率跌得愈少, 就會輸給大盤. 當然在台灣買, 還有匯率變動問題. 因此買長債這個送分題只有邏輯正確, 川普都想讓你賺, 但這次時間之神還沒同意.