3-3乗加重関数（tri-cubic weighting function）

tri_cubic <- function(x) {
  ifelse(abs(x) <= 1, (1 - abs(x)^3)^3, 0)
}
x_vals <- seq(-10, 10, length.out = 47)
x_vals <- seq(-10, 10, length.out = 100)
y_vals <- sapply(x_vals, tri_cubic)         

plot(x_vals, y_vals, type = "l", col = "skyblue", lwd = 2,
     xlab = "x", ylab = "w(x)", main = "Tri-Cubic Weighting Function")
abline(h = 0, v = 0, col = "gray", lty = 100)  # 補助線を追加

12月15日から予測せよ。

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

## 要求されたパッケージ quantmod をロード中です

## Warning: パッケージ 'quantmod' はバージョン 4.4.2 の R の下で造られました

## 要求されたパッケージ xts をロード中です

## Warning: パッケージ 'xts' はバージョン 4.4.2 の R の下で造られました

## 要求されたパッケージ zoo をロード中です

## Warning: パッケージ 'zoo' はバージョン 4.4.2 の R の下で造られました

## 
## 次のパッケージを付け加えます: 'zoo'

## 以下のオブジェクトは 'package:base' からマスクされています:
## 
##     as.Date, as.Date.numeric

## 要求されたパッケージ TTR をロード中です

## Warning: パッケージ 'TTR' はバージョン 4.4.2 の R の下で造られました

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

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)
}

# 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

##              i    y
## 2022-01-04   1 5482
## 2022-01-05   2 5420
## 2022-01-06   3 5372
## 2022-01-07   4 5489
## 2022-01-11   5 5360
## 2022-01-12   6 5683
## 2022-01-13   7 5566
## 2022-01-14   8 5497
## 2022-01-17   9 5530
## 2022-01-18  10 5508
## 2022-01-19  11 5412
## 2022-01-20  12 5527
## 2022-01-21  13 5487
## 2022-01-24  14 5355
## 2022-01-25  15 5069
## 2022-01-26  16 5156
## 2022-01-27  17 4692
## 2022-01-28  18 4795
## 2022-01-31  19 5011
## 2022-02-01  20 5067
## 2022-02-02  21 5201
## 2022-02-03  22 5195
## 2022-02-04  23 5214
## 2022-02-07  24 5350
## 2022-02-08  25 5302
## 2022-02-09  26 5612
## 2022-02-10  27 5483
## 2022-02-14  28 5268
## 2022-02-15  29 5158
## 2022-02-16  30 5236
## 2022-02-17  31 5117
## 2022-02-18  32 5184
## 2022-02-21  33 5160
## 2022-02-22  34 5154
## 2022-02-24  35 4802
## 2022-02-25  36 5069
## 2022-02-28  37 5125
## 2022-03-01  38 5232
## 2022-03-02  39 5266
## 2022-03-03  40 5212
## 2022-03-04  41 4963
## 2022-03-07  42 4707
## 2022-03-08  43 4474
## 2022-03-09  44 4636
## 2022-03-10  45 4780
## 2022-03-11  46 4483
## 2022-03-14  47 4450
## 2022-03-15  48 4265
## 2022-03-16  49 4519
## 2022-03-17  50 4785
## 2022-03-18  51 4961
## 2022-03-22  52 5056
## 2022-03-23  53 5421
## 2022-03-24  54 5498
## 2022-03-25  55 5402
## 2022-03-28  56 5402
## 2022-03-29  57 5504
## 2022-03-30  58 5632
## 2022-03-31  59 5559
## 2022-04-01  60 5578
## 2022-04-04  61 5785
## 2022-04-05  62 5936
## 2022-04-06  63 5769
## 2022-04-07  64 5654
## 2022-04-08  65 5683
## 2022-04-11  66 5530
## 2022-04-12  67 5482
## 2022-04-13  68 5592
## 2022-04-14  69 5762
## 2022-04-15  70 5692
## 2022-04-18  71 5653
## 2022-04-19  72 5547
## 2022-04-20  73 5627
## 2022-04-21  74 5579
## 2022-04-22  75 5411
## 2022-04-25  76 4989
## 2022-04-26  77 5195
## 2022-04-27  78 5182
## 2022-04-28  79 5269
## 2022-05-02  80 5290
## 2022-05-06  81 5170
## 2022-05-09  82 4989
## 2022-05-10  83 4900
## 2022-05-11  84 4883
## 2022-05-12  85 4491
## 2022-05-13  86 5040
## 2022-05-16  87 5110
## 2022-05-17  88 5122
## 2022-05-18  89 5137
## 2022-05-19  90 5055
## 2022-05-20  91 5232
## 2022-05-23  92 5276
## 2022-05-24  93 5191
## 2022-05-25  94 5102
## 2022-05-26  95 5164
## 2022-05-27  96 5341
## 2022-05-30  97 5358
## 2022-05-31  98 5379
## 2022-06-01  99 5300
## 2022-06-02 100 5302
## 2022-06-03 101 5418
## 2022-06-06 102 5391
## 2022-06-07 103 5430
## 2022-06-08 104 5563
## 2022-06-09 105 5659
## 2022-06-10 106 5545
## 2022-06-13 107 5165
## 2022-06-14 108 5033
## 2022-06-15 109 5028
## 2022-06-16 110 5000
## 2022-06-17 111 4788
## 2022-06-20 112 4933
## 2022-06-21 113 5075
## 2022-06-22 114 5025
## 2022-06-23 115 5103
## 2022-06-24 116 5224
## 2022-06-27 117 5418
## 2022-06-28 118 5414
## 2022-06-29 119 5325
## 2022-06-30 120 5235
## 2022-07-01 121 5131
## 2022-07-04 122 5283
## 2022-07-05 123 5378
## 2022-07-06 124 5343
## 2022-07-07 125 5388
## 2022-07-08 126 5400
## 2022-07-11 127 5446
## 2022-07-12 128 5213
## 2022-07-13 129 5338
## 2022-07-14 130 5385
## 2022-07-15 131 5296
## 2022-07-19 132 5371
## 2022-07-20 133 5500
## 2022-07-21 134 5507
## 2022-07-22 135 5482
## 2022-07-25 136 5505
## 2022-07-26 137 5681
## 2022-07-27 138 5688
## 2022-07-28 139 5656
## 2022-07-29 140 5605
## 2022-08-01 141 5504
## 2022-08-02 142 5450
## 2022-08-03 143 5496
## 2022-08-04 144 5640
## 2022-08-05 145 5653
## 2022-08-08 146 5695
## 2022-08-09 147 5295
## 2022-08-10 148 5315
## 2022-08-12 149 5610
## 2022-08-15 150 5900
## 2022-08-16 151 5749
## 2022-08-17 152 5802
## 2022-08-18 153 5815
## 2022-08-19 154 5764
## 2022-08-22 155 5739
## 2022-08-23 156 5600
## 2022-08-24 157 5616
## 2022-08-25 158 5638
## 2022-08-26 159 5705
## 2022-08-29 160 5496
## 2022-08-30 161 5588
## 2022-08-31 162 5562
## 2022-09-01 163 5510
## 2022-09-02 164 5493
## 2022-09-05 165 5497
## 2022-09-06 166 5479
## 2022-09-07 167 5368
## 2022-09-08 168 5499
## 2022-09-09 169 5500
## 2022-09-12 170 5637
## 2022-09-13 171 5654
## 2022-09-14 172 5405
## 2022-09-15 173 5438
## 2022-09-16 174 5458
## 2022-09-20 175 5517
## 2022-09-21 176 5453
## 2022-09-22 177 5343
## 2022-09-26 178 5066
## 2022-09-27 179 5030
## 2022-09-28 180 4938
## 2022-09-29 181 5031
## 2022-09-30 182 4900
## 2022-10-03 183 4981
## 2022-10-04 184 5235
## 2022-10-05 185 5355
## 2022-10-06 186 5494
## 2022-10-07 187 5504
## 2022-10-11 188 5451
## 2022-10-12 189 5500
## 2022-10-13 190 5469
## 2022-10-14 191 5650
## 2022-10-17 192 5549
## 2022-10-18 193 5537
## 2022-10-19 194 5743
## 2022-10-20 195 5780
## 2022-10-21 196 5758
## 2022-10-24 197 5733
## 2022-10-25 198 5949
## 2022-10-26 199 6031
## 2022-10-27 200 6087
## 2022-10-28 201 6026
## 2022-10-31 202 6400
## 2022-11-01 203 6574
## 2022-11-02 204 6705
## 2022-11-04 205 6557
## 2022-11-07 206 6601
## 2022-11-08 207 6929
## 2022-11-09 208 7019
## 2022-11-10 209 6833
## 2022-11-11 210 6953
## 2022-11-14 211 6068
## 2022-11-15 212 6147
## 2022-11-16 213 6328
## 2022-11-17 214 6319
## 2022-11-18 215 6075
## 2022-11-21 216 6077
## 2022-11-22 217 6068
## 2022-11-24 218 6064
## 2022-11-25 219 6079
## 2022-11-28 220 6042
## 2022-11-29 221 5959
## 2022-11-30 222 5952
## 2022-12-01 223 6047
## 2022-12-02 224 6055
## 2022-12-05 225 6068
## 2022-12-06 226 5917
## 2022-12-07 227 5967
## 2022-12-08 228 6095
## 2022-12-09 229 6168
## 2022-12-12 230 6205
## 2022-12-13 231 6213
## 2022-12-14 232 6243
## 2022-12-15 233 6291
## 2022-12-16 234 6052
## 2022-12-19 235 6018
## 2022-12-20 236 5726
## 2022-12-21 237 5800
## 2022-12-22 238 5786
## 2022-12-23 239 5757
## 2022-12-26 240 5778
## 2022-12-27 241 5796
## 2022-12-28 242 5711
## 2022-12-29 243 5618
## 2022-12-30 244 5644

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

xt.tr <- xt["2022-01-01/2022-12-14", is.close] # 11月迄の株価（終値）
xt.te <- xt["2022-12-15/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))

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％が含まれる。

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シグマ区間"))

予測値が上向きになっているのに対し、実際の株価は下落してしまっている