1 局所回帰（local regression）

1.1 LOESS
（locally estimated scattered smoothing）

代表的な局所回帰であるLOESS（ロエス）について説明する。

1.2 平滑化窓幅

\(平滑化窓幅 = \frac{近傍点数k}{標本サイズn}\)が設定すべきパラメータである。 平滑化窓幅（span）を決めると標本サイズ\(n\)から近傍点数\(k\)が決まる。 平滑化したい位置（時刻）\(x_i\)から\(x\)軸上の距離的（時間的） に最も近い\(k\)個の標本点が選ばれる。

1.3 加重関数

3-3乗加重関数（tri-cubic weighting function）とは， 近傍ほど加重が大きい，つまり重要とみなす関数である。 これを使用し回帰係数を推定するための加重を計算する。

\[ f(d)= \begin{cases} \left(1-|d|^3\right)^3\quad |d|<1\\ \hspace{8mm}0\hspace{16mm} |d|\ge 1 \end{cases} \]

ここで，\(d = \frac{近傍点までの距離}{近傍点までの距離の中で最大の距離},\quad -1 < d < 1\)とする規格化距離。

1.4 多項式回帰モデル
（polynomial regression model）

LOESSでは，\(k\)個の近傍点データに対して， 1次や2次の比較的低次の多項式回帰モデルをフィッティングさせる。 \(j\)を近傍点のインデックス，\(\beta_0\)，\(\beta_1\)，…を（偏）回帰係数をとすると， 次のようなモデルとなる。

  • 1次多項式回帰（単回帰）モデル

\[Y_j = \beta_0 + \beta_1 x_j + \epsilon_j, \quad \epsilon_j \sim \mathbf{N}(0, \sigma^2)\]

  • 2次多項式回帰モデル（デフォルト）

\[Y_j = \beta_0 + \beta_1 x_j + \beta_2 x_j^2 + \epsilon_j,\quad \epsilon_j \sim \mathbf{N}(0, \sigma^2)\]

\(k\)個の近傍点データを用いて， 次の目的関数\(L\)を最小とする回帰係数を求める。 位置（時刻）\(x_i\)が決まれば， 近傍点への距離\(d_j\)が得られ，加重\(w_j=f(d_j)\)が決まる。 \[ 目的関数：　L = \sum_j w_j(y_j - \hat{y_j})^2 \] ここで，\(y_j\)，\(\hat{y_j}\)は，それぞれ近傍点値，回帰モデルによる予測値である。

回帰係数の最小２乗推定値を\(b_0, b_1, b_2\)とすると， 平滑値\(s_i\)は次のようになる。

\[s_i = b_0 + b_1 x_i\] \[s_i = b_0 + b_1 x_i + b_2 x_i^2\]

2 実データを使った分析例

2.1 株価のトレンド分析

ヤフーファイナンス(yahooj)より携帯電話会社の株価を取得しトレンド分析する。

何回もサイトにアクセスしサーバーに負荷をかけないように， 1度ダウンロードしたファイルを使う。
新規に取得したい場合は保存ファイルを削除すること。

if (!require(quantmod)) install.packages("quantmod")

symbol <- "9984" # 証券コード（ソフトバンク）
#symbol <- "4755" # 証券コード（楽天）
#symbol <- "9433" # 証券コード（AU）

FILE <- paste0(symbol, ".csv") # 保存用CSVファイル
MAIN <- paste0("株価（証券コード：", symbol, "）") # 主タイトルラベル

if ( !file.exists(FILE) )
{ 
  # CSVファイルが存在しない場合，サイトからデータをダウンロードする。
  xt <- getSymbols(symbol,
                   src  = "yahooj",
                   from = "2022-01-01",
                   to   = "2022-12-31",
                   auto.assign = F)

   # 株価データをCSVファイルとして保存する。
   write.zoo(xt, file = FILE, sep = ",", quote = F)
   
} else
{
  # CSVファイルが存在する場合，そこから株価データを取得
  xt <- getSymbols(symbol, src = "csv", auto.assign = F)
}

2.1.1 インタラクティブグラフ

if (!require(dygraphs)) install.packages("dygraphs")

# カラム名からティッカー名である接頭語（.の前まで）を除去
colnames(xt) <- substring(colnames(xt), regexpr("\\.[^\\.]*$", colnames(xt)) + 1)

xt |> HLC() |>
  dygraph(main = "ソフトバンク（9984）") |>
  dyAxis("y", label = "Yen") |>
  dySeries(c("Low", "Close", "High")) |>
  dyRangeSelector()
xt |> tail(n = 30) |> OHLC() |>
  dygraph(main = "ソフトバンク（9984）") |>
  dyAxis("y", label = "Yen") |>
  dyCandlestick()
# Closeという文字列を含むカラムを探す。
# CSVから読み込むと何故かカラム名からYJの文字が消えるためこのようにしている。
is.close <- grepl("Close", colnames(xt)) # 終値（Close）カラムか調べる。
y <- as.vector(xt[, is.close])           # 終値データをベクトルとして取得

x <- as.POSIXct(index(xt)) # タイムスタンプ
i <- seq_along(x)          # タイムスタンプのインデックス
n <- length(x)             # 標本サイズ（sample size）
d <- data.frame(i, y)
rownames(d) <- x
d

2.2 LOESSによるフィッティング

# 多項式の次数 (degree)はデフォルト2
# 窓幅 (span) = 近傍点数/標本サイズはデフォルト0.75
fit10 <- loess(y ~ i, data = d, span = 0.1)
fit30 <- loess(y ~ i, data = d, span = 0.3)

2.3 グラフ

# 図枠 (diagram frame)
matplot(x, y, type = "n", main = MAIN, xlab = "日付", ylab = "円")

# 格子線 (grid lines)
x.g <- seq(as.POSIXct("2020-01-01"),
           as.POSIXct("2030-01-01"), by = 'month')

abline(h = seq(0, 9000, 500), v = x.g, lty = 2, lwd = 0.5, col = "lightgray")

# 時系列線 (time series lines)
matlines(x, y,            col = 1, lwd = 2) # 原系列 (raw series)
matlines(x, fit10$fitted, col = 2, lwd = 2) # 窓幅 10％
matlines(x, fit30$fitted, col = 3, lwd = 2) # 窓幅 30％

# 凡例 (legend)
legend("topleft", col = 1:3, lty = 1,
       legend = c("株価", "LOESS（窓幅 10%）", "LOESS（窓幅 30%）"))

3 予測

3.1 データ

訓練データとして1〜11月株価（終値）を使用し10日間予測する。 正解データとして12月株価（終値）を使用する。

xt.tr <- xt["2022-01-01/2022-11-30", is.close] # 11月迄の株価（終値）
xt.te <- xt["2022-12-01/2022-12-31", is.close] # 12月株価（終値）
xt.al <- xt["2022-01-01/2022-12-31", is.close] # 2022年株価（終値）

x.tr <- as.POSIXct(index(xt.tr))
x.te <- as.POSIXct(index(xt.te))
x.al <- as.POSIXct(index(xt.al))

n.tr <- length(x.tr)
n.te <- length(x.te)
n.al <- length(x.al)

d.tr <- data.frame(x = x.tr, i = 1:n.tr, y = as.vector(xt.tr))
d.te <- data.frame(x = x.te, i = 1:n.te, y = as.vector(xt.te))
d.al <- data.frame(x = x.al, i = 1:n.al, y = as.vector(xt.al))

3.2 LOESSによるフィッティング

多項式の次数 (degree)はデフォルト2
平滑化窓幅 (span) = 近傍点数/標本サイズはデフォルト0.75
企業ごとに，その株価の変化の激しさに応じてこれらのパラメータを調整する必要がある。

SPAN <- 0.1

fit <- loess(y ~ i, data = d.tr, span = SPAN,
             control = loess.control(surface = "direct"))
N.DAYS.AHEAD <- 10 # 予測日数 (forecasting horizon in days)

pred <- predict(fit, se = T,
                newdata = data.frame(i = (n.tr+1):(n.tr+N.DAYS.AHEAD)))

yhat <- pred$fit           # 予測値 (predicted values)
upp  <- yhat + pred$se.fit # 1シグマ区間（上限）
low  <- yhat - pred$se.fit # 1シグマ区間（下限）
# 正規分布を仮定すると閉区間[low, upp]にデータの68％が含まれる。

3.3 グラフ

SEC.30DAYS <- 60*60*24*30 # 30日の秒数 

# 図枠 (diagram frame)
matplot(x = x.al,
        y = d.al$y,
        type = "n", # プロットしない。 
        xaxt = "n", # x軸を描かない。
        yaxt = "n", # y軸を描かない。
        xlim = c(x.te[1] - 2*SEC.30DAYS, x.te[n.te]), # 画面表示期間
        ylim = c(0.8*min(d.al$y, low), 1.2*max(d.al$y, upp)),
        main = MAIN,
        xlab = paste(unique(format(x.al, "%Y")), "年"),
        ylab = "円")

y.g <- seq(0, 10000, 500)
axis(side = 1, tick = F, at = x.al,
     labels = format(x.al,  "%m/%d(%a)")) # x軸
axis(side = 2, tick = T, at = y.g,
     labels = format(y.g, big.mark = ",")) # y軸

# 格子線 (grid lines)
x.g <- seq(as.POSIXct("2020-01-01"),
            as.POSIXct("2030-01-01"), by = "weeks")

abline(h = seq(0, 9000, 500), v = x.g, lty = 2, lwd = 0.5, col = "lightgray")

# 時系列線 (time series lines)
COL <- c("gray", "black", "blue", "darkgreen")
PCH <- c(1, 2, 16, NA)
LTY <- c(rep(1, 3), 2)

# プロット
matlines(x.al, d.al$y,
         type = "o", pch = PCH[1], col = COL[1], lty = LTY[1])
matlines(x.tr, fit$fitted,
         type = "o", pch = PCH[2], col = COL[2], lty = LTY[2])
matlines(x.te[1:N.DAYS.AHEAD], yhat,
         type = "o", pch = PCH[3], col = COL[3], lty = LTY[3])
matlines(x.te[1:N.DAYS.AHEAD], upp,
         type = "l", pch = PCH[4], col = COL[4], lty = LTY[4])
matlines(x.te[1:N.DAYS.AHEAD], low,
         type = "l", pch = PCH[4], col = COL[4], lty = LTY[4])

# 予測開始線 (start line for forecasting)
abline(v = x.tr[n.tr] + 60*60*12, col = 3, lty = 2) 

# 凡例 (legend)
legend("topleft", pch = PCH, col = COL, lty = LTY,
       legend = c("株価（終値）",
                  sprintf("平滑値（窓幅%3.0d%%）", SPAN*100),
                  "予測値",
                  "1シグマ区間"))

4 演習課題