機械学習

演習課題

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

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

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

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

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

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

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

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

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

FILE <- paste0(symbol, ".csv") 
MAIN <- paste0("株価（証券コード：", symbol, "）") 

if ( !file.exists(FILE) )
{ 
  
  xt <- getSymbols(symbol,
                   src  = "yahooj",
                   from = "2022-01-01",
                   to   = "2022-12-31",
                   auto.assign = F)

   write.zoo(xt, file = FILE, sep = ",", quote = F)
   
} else
{
  xt <- getSymbols(symbol, src = "csv", auto.assign = F)
}

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

##  要求されたパッケージ 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()

is.close <- grepl("Close", colnames(xt))

y <- as.vector(xt[, is.close])           

x <- as.POSIXct(index(xt)) 
i <- seq_along(x)          
n <- length(x)             
d <- data.frame(i, y)
rownames(d) <- x
d

fit10 <- loess(y ~ i, data = d, span = 0.1)
fit30 <- loess(y ~ i, data = d, span = 0.3)

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) 
matlines(x, fit10$fitted, col = 2, lwd = 2) 
matlines(x, fit30$fitted, col = 3, lwd = 2) 

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

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

SPAN <- 0.1

fit <- loess(y ~ i, data = d.tr, span = SPAN,
             control = loess.control(surface = "direct"))

N.DAYS.AHEAD <- 10 

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

yhat <- pred$fit           
upp  <- yhat + pred$se.fit 
low  <- yhat - pred$se.fit

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

abline(v = x.tr[n.tr] + 60*60*12, col = 3, lty = 2) 

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

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

xt.tr <- xt["2022-01-01/2022-12-14", is.close]
xt.te <- xt["2022-12-15/2022-12-31", is.close] 
xt.al <- xt["2022-01-01/2022-12-31", is.close]

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 

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

yhat <- pred$fit           
upp  <- yhat + pred$se.fit 
low  <- yhat - pred$se.fit

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

abline(v = x.tr[n.tr] + 50*50*12, col = 3, lty = 2) 

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