f <- function(d) ifelse(abs(d) < 1, (1 - abs(d)^3)^3, 0)
d <- seq(-2, 2, 0.1)
matplot(x = d, y = f(d), type = 'l', col = 4, main = '3-3重加重関数(tri-cubic weighting function)')
grid()

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)
}
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
# 多項式の次数 (degree)はデフォルト2
# 窓幅 (span) = 近傍点数/標本サイズはデフォルト0.75
fit10 <- loess(y ~ i, data = d, span = 0.1)
fit30 <- loess(y ~ i, data = d, span = 0.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%)"))

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

予測値は上昇傾向になっているのに対し、実際の株価は下落している。教訓として、市場の動きが必ずしも予測通りになるとは限らない。