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

1 移動平均

次のロボット販売のデータを使って，移動平均を求めトレンドを分析する。

1.1 データ

ロボット販売台数

year  <- 2011:2017
sales <- c(1, 3, 2, 7, 3, 8, 7)
n <- length(year)
d <- data.frame(year, sales)
d

1.2 グラフ

plot_ly(x = d$year, y = d$sales, type = "scatter", mode = "markers+lines") |>
  layout(title = "ロボット販売",
         xaxis = list(title = "会計年"),
         yaxis = list(title = "台数"))

講義資料（PDF）

cairo_pdf('robot_sales.pdf')
matplot(x = year, y = sales, type = 'o', pch = 1, col = 4,
        main = 'ロボット販売',
        xlab = '会計年',
        ylab = '台数')
grid()
dev.off()
## png 
##   2

2 移動平均の計算

2.1 項数が奇数の場合

項数3の移動平均\(s_i\)を求める。 なお，移動する計算範囲の両端で原系列が存在しないとき， \(s_i\)は計算できないため\(NA\) (not available)となる。

\[ s_i = (y_{i-1} + y_i + y_{i+1}) / 3 \]

d$ma3 <- NA # 値がNAの移動平均のカラム（ma3）を追加

for (i in 2:(n-1))
{
  d$ma3[i] <- (sales[i-1] + sales[i] + sales[i+1]) / 3
}

d$ma3
## [1] NA  2  4  4  6  6 NA

2.2 Rパッケージを利用する場合

# filter関数を使用した（中心化）移動平均
ma3 <- stats::filter(sales, rep(1/3, 3))
ma3
## Time Series:
## Start = 1 
## End = 7 
## Frequency = 1 
## [1] NA  2  4  4  6  6 NA
# quantmodパッケージのrollmean関数を使用する場合
ma3 <- rollmean(sales, 3, na.pad = T)
ma3
## [1] NA  2  4  4  6  6 NA

2.3 項数が偶数の場合

項数4の移動平均\(\hat{y_i}\) を求める。 1項多い5項にし両端の値の重みを0.5にして全部で4項相当にして計算する。
結局，奇数項数にして計算するので初めから奇数項数で設定したほうがよい。

\[ s_i = (0.5\times y_{i-2} + y_{i-1} + y_i + y_{i+1} + 0.5\times y_{i+2}) / 4 \]

d$ma4 <- NA # 値がNAの移動平均のカラム（ma4）を追加

for (i in 3:(n-2)) 
{
  d$ma4[i] <-
    (0.5 * sales[i-2] + sales[i-1] + sales[i] + sales[i+1] + 0.5 * sales[i+2]) / 4
}

2.4 Rパッケージを利用する場合

# filter関数を使用する場合
  ma4 <- stats::filter(sales, c(0.5/4, 1/4, 1/4, 1/4, 0.5/4))
  ma4
## Time Series:
## Start = 1 
## End = 7 
## Frequency = 1 
## [1]    NA    NA 3.500 4.375 5.625    NA    NA
# quantmodパッケージのrollmean関数を使用する場合
# 計算方法が異なり単純な4項の平均であるため中心が1つずれるので注意
  ma4 <- rollmean(sales, k = 4, na.pad = T)
  ma4
## [1]   NA 3.25 3.75 5.00 6.25   NA   NA

2.5 移動平均まとめ

中心化移動平均（ma3：次数3，ma4：次数4）

options(digits = 2)
d

2.6 グラフ

plot_ly(type = "scatter", mode = "markers+lines") |>
  add_trace(x = d$year, y = d$sales, name = "台数") |>
  add_trace(x = d$year, y = d$ma3,   name = "3項移動平均") |>
  add_trace(x = d$year, y = d$ma4,   name = "4項移動平均") |>
  layout(title = "ロボット販売",
         xaxis = list(title = "会計年"),
         yaxis = list(title = "台数"))

講義資料（PDF）

cairo_pdf('robot_sales_moving_average.pdf')

matplot(x = year, y = sales, type = 'o', pch = 1,
        main = 'ロボット販売', xlab = '会計年', ylab = '台数')
grid()

matlines(x = year, y = d$ma3, type = 'o', pch = 2, col = 2)
matlines(x = year, y = d$ma4, type = 'o', pch = 4, col = 4)

legend('topleft', col = c(1, 2, 4), pch = c(1, 2, 4), lty = 1,
       legend = c('台数', '3項移動平均', '4項移動平均'))

dev.off()
## png 
##   2

3 実データを用いたトレンド分析例

3.1 新型コロナウィルス感染症の陽性者数のトレンド分析

コロナ陽性者数のトレンドを31日移動平均で行う。

【厚生労働省】新規陽性者数の推移（日別）

d0 <- read.csv(file = 'newly_confirmed_cases_daily.csv')
d0
x <- as.POSIXct(d0$Date, format = '%Y/%m/%d') # 「/」で区切られていることに注意
y <- d0$ALL

yhat <- rollmean(y, 31, na.pad = T)

plot_ly(type = "scatter", mode = "lines") |>
  add_trace(x = x, y = y,    name = "原系列") |>
  add_trace(x = x, y = yhat, name = "31日移動平均") |>
  config(locale = "ja") |> # x軸ラベルの月を日本語表示
  layout(title = "新規陽性者数の推移（日別）",
         xaxis = list(title = "年"),
         yaxis = list(title = "陽性者数"))

講義資料（PDF）

cairo_pdf('corona.pdf')

matplot(x = x, y = y, type = 'l', col = 3, 
        main = '新規陽性者数の推移（日別）', xlab = '年', ylab = '陽性者数')
matlines(x = x, y = yhat, col = 1)
grid()
legend('topleft', col = c(3, 1), lty = 1, legend = c('原系列', '31日移動平均'))

dev.off()
## png 
##   2

3.2 株価のトレンド分析

株価のテクニカル分析で使用する移動平均は， 当日のデータは得られない 状況下で当日の予測値が必要となるため次の定義となる。 ただし，この移動平均は時間遅れ（time lag）が生じるので要注意。

（３日移動平均） \[ \begin{align} s_i &= (3日前終値 + 2日前終値 + 1日前終値) / 3\\ &= (y_{i-3} + y_{i-2} +y_{i-1}) / 3 \end{align} \] ソフトバンクの証券コード：9984
ヤフーファイナンス(yahooj)よりデータを取得

z <- getSymbols('9984', src  = 'yahooj',
                        from = '2022-01-01',
                        to   = '2022-12-31',
                        auto.assign = F)

# データの保存
write.zoo(z, file = 'softbank.csv', sep = ',', quote = F)

x <- as.POSIXct(index(z), format = '%Y-%m-%d')
y <- z[, 'YJ9984.Close']

# 1日分データをずらすことで当日の日付の株価が１日前のものになる。
ylag1 <- lag(y, k = 1)

# align = 'right'で右寄せ移動平均となる。
yhat5   <- rollmean(ylag1,  5, fill = NA, align = 'right')
yhat25  <- rollmean(ylag1, 25, fill = NA, align = 'right')
yhat25c <- rollmean(y, 25, fill = NA) # 参考（通常の移動平均：中心化移動平均）

matplot(x, y, type = 'l', col = 3,
main = 'ソフトバンク（9984）', xlab = '日', ylab = '円')
matlines(x, yhat5,   col = 2)
matlines(x, yhat25,  col = 4)
matlines(x, yhat25c, col = 1)
grid()
legend('topleft', col = c(3, 2, 4, 1), lty = 1,
       legend = c('株価',
                  '5日移動平均',
                  '25日移動平均',
                  '25日中心化移動平均（参考）'))

# （参考）quantmodパッケージのグラフ機能
chartSeries(z, TA = c(addBBands(), addMACD()), type='bar',
            subset='2022-06-01/2022-12-31')

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

y値がxtsオブジェクトになっており， plot_ly関数はこの形式のオブジェクトを受け付けないので paste関数で文字列に変換する必要がある。

plot_ly(type = "scatter", mode = "lines") |>
  add_trace(x = x, y = paste(y),       name = "株価") |>
  add_trace(x = x, y = paste(yhat5),   name = "5日移動平均") |>
  add_trace(x = x, y = paste(yhat25),  name = "25日移動平均") |>
  add_trace(x = x, y = paste(yhat25c), name = "25日中心化移動平均（参考）") |>
  config(locale = "ja") |> # x軸ラベルの月を日本語表示
  layout(title = "ソフトバンク（9984）",
         xaxis = list(title = "日"),
         yaxis = list(title = "円"))

3.3 局所多項式回帰（LOESS: locally estimated scattered smoothing）

時間遅れが無いことに注目

# LOESS（平滑化窓幅: span）
x <- seq_along(y)
fit10 <- loess(y ~ x, data = data.frame(x, y), span = 0.1) # 10%
fit30 <- loess(y ~ x, data = data.frame(x, y), span = 0.3) # 30%

matplot(x, y = y, type = "l", col = 3,
main = "ソフトバンク（9984）", xlab = "日", ylab = "円")
matlines(x, fit10$fitted, col = 1, lwd = 2)
matlines(x, fit30$fitted, col = 2, lwd = 2)
grid()
legend("topleft", col = c(3, 1, 2), lty = 1,
       legend = c("株価",
                  "LOESS（平滑化窓幅10%）",
                  "LOESS（平滑化窓幅30%）"))

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

y値がxtsオブジェクトになっており， plot_ly関数はこの形式のオブジェクトを受け付けないので paste関数で文字列に変換する必要がある。

plot_ly(type = "scatter", mode = "lines") |>
  add_trace(x = x, y = paste(y),            name = "株価") |>
  add_trace(x = x, y = paste(fit10$fitted), name = "LOESS（平滑化窓幅10%）") |>
  add_trace(x = x, y = paste(fit30$fitted), name = "LOESS（平滑化窓幅30%）") |>
  layout(title = "ソフトバンク（9984）", 
         xaxis = list(title = "日"),
         yaxis = list(title = "円"))

4 演習課題

  1. 【厚生労働省】PCR検査実施人数 のトレンド分析を行え

  2. ZOZO（証券コード：3092）の株価トレンド分析を行え