機械学習
演習課題
2023-11-21
- このページは公開されるので個人情報などは記載しないこと。
- 課題が完成したらknitしPublishする。
課題名
rm(list = ls()) # 全オブジェクト削除
# ここにRコードを記入する。
#例）
plot(1:9)
#時系列分析(局所回帰)
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()
library(quantmod)## Warning: パッケージ 'quantmod' はバージョン 4.3.2 の R の下で造られました## 要求されたパッケージ xts をロード中です## Warning: パッケージ 'xts' はバージョン 4.3.2 の R の下で造られました## 要求されたパッケージ zoo をロード中です## Warning: パッケージ 'zoo' はバージョン 4.3.2 の R の下で造られました##
## 次のパッケージを付け加えます: 'zoo'## 以下のオブジェクトは 'package:base' からマスクされています:
##
## as.Date, as.Date.numeric## 要求されたパッケージ TTR をロード中です## Warning: パッケージ 'TTR' はバージョン 4.3.2 の R の下で造られました## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
Symbols <- '9984' # 証券コード（ソフトバンク）
#Symbols <- '4755' # 証券コード（楽天）
#Symbols <- '9433' # 証券コード（AU）
FILE <- paste0(Symbols, '.csv') # 保存用CSVファイル
MAIN <- paste0('株価（証券コード：', Symbols, '）') # 主タイトルラベル
if ( !file.exists(FILE) )
{
# CSVファイルが存在しない場合，サイトからデータをダウンロードする。
xt <- getSymbols(Symbols,
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(Symbols, src = 'csv', auto.assign = F)
}
# Closeという文字列を含むカラムを探す。
# CSVから読み込むと何故かカラム名からYJの文字が消えるためこのようにしている。
is.close <- grepl('Close', colnames(xt)) # 終値（Close）カラムか調べる。
y <- as.vector(xt[, is.close]) # 終値データをベクトルとして取得
px <- as.POSIXct(index(xt)) # タイムスタンプ
i <- seq_along(px) # タイムスタンプのインデックス
n <- length(px) # 標本サイズ（sample size）
d <- data.frame(i, y)
rownames(d) <- px
library(DT)## Warning: パッケージ 'DT' はバージョン 4.3.2 の R の下で造られましたdatatable(d, caption = MAIN, colnames = c('連番', '終値（円）'))# 多項式の次数 (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(px, y, type = 'n',
main = MAIN,
ylab = '円')
# 格子線 (grid lines)
px.g <- seq(as.POSIXct('2020-01-01'),
as.POSIXct('2030-01-01'), by = 'month')
abline(h = seq(0, 9000, 500), v = px.g, lty = 2, lwd = 0.5, col = 'lightgray')
# 時系列線 (time series lines)
matlines(px, y, col = 1, lwd = 2) # 原系列 (raw series)
matlines(px, fit10$fitted, col = 2, lwd = 2) # 窓幅 10％
matlines(px, fit30$fitted, col = 3, lwd = 2) # 窓幅 30％
# 凡例 (legend)
legend('topright', 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年株価（終値）
px.tr <- as.POSIXct(index(xt.tr))
px.te <- as.POSIXct(index(xt.te))
px.al <- as.POSIXct(index(xt.al))
n.tr <- length(px.tr)
n.te <- length(px.te)
n.al <- length(px.al)
d.tr <- data.frame(px = px.tr, i = 1:n.tr, y = as.vector(xt.tr))
d.te <- data.frame(px = px.te, i = 1:n.te, y = as.vector(xt.te))
d.al <- data.frame(px = px.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 = px.al,
y = d.al$y,
type = 'n', # プロットしない。
xaxt = 'n', # x軸を描かない。
yaxt = 'n', # y軸を描かない。
xlim = c(px.te[1] - 2*SEC.30DAYS, px.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(px.al, '%Y')), '年'),
ylab = '円')
y.g <- seq(0, 10000, 500)
axis(side = 1, tick = F, at = px.al, labels = format(px.al, '%m/%d(%a)')) # x軸
axis(side = 2, tick = T, at = y.g, labels = format(y.g, big.mark = ',')) # y軸
# 格子線 (grid lines)
px.g <- seq(as.POSIXct('2020-01-01'),
as.POSIXct('2030-01-01'), by = 'weeks')
abline(h = seq(0, 9000, 500), v = px.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(px.al, d.al$y, type = 'o', pch = PCH[1], col = COL[1], lty = LTY[1])
matlines(px.tr, fit$fitted, type = 'o', pch = PCH[2], col = COL[2], lty = LTY[2])
matlines(px.te[1:N.DAYS.AHEAD], yhat, type = 'o', pch = PCH[3], col = COL[3], lty = LTY[3])
matlines(px.te[1:N.DAYS.AHEAD], upp, type = 'l', pch = PCH[4], col = COL[4], lty = LTY[4])
matlines(px.te[1:N.DAYS.AHEAD], low, type = 'l', pch = PCH[4], col = COL[4], lty = LTY[4])
# 予測開始線 (start line for forecasting)
abline(v = px.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シグマ区間'))
#演習
xt.tr <- xt['2022-01-01/2022-12-15', is.close] # 11月迄の株価（終値）
xt.te <- xt['2022-12-01/2022-12-15', is.close] # 12月株価（終値）
xt.al <- xt['2022-01-01/2022-12-31', is.close] # 2022年株価（終値）
px.tr <- as.POSIXct(index(xt.tr))
px.te <- as.POSIXct(index(xt.te))
px.al <- as.POSIXct(index(xt.al))
n.tr <- length(px.tr)
n.te <- length(px.te)
n.al <- length(px.al)
d.tr <- data.frame(px = px.tr, i = 1:n.tr, y = as.vector(xt.tr))
d.te <- data.frame(px = px.te, i = 1:n.te, y = as.vector(xt.te))
d.al <- data.frame(px = px.al, i = 1:n.al, y = as.vector(xt.al))