【forecast】Rでファイナンスデータの時系列分析【ARIMA】

| 1件のコメント

計量時系列分析について調べたこと, 主に AR, ARMA, ARIMA モデルについてメモります。

ファイナンスデータの準備

今回は Yahoo!ファイナンス から日産自動車[7201]のデータを取得しました。

csvの文字コードを Shift_JIS から utf-8 に変更しておきます。OS Xであれば nkfコマンド。

$ brew install nkf
$ nkf -w --overwrite ./data/7201.T.csv

まずは, {quantmod} でロウソク足を図示してみます。

require(quantmod)

# convert char-code "nkf -w --overwrite filename"
df <- read.csv("data/7201.T.csv", header = T)
names(df) <- c("Date", "Open", "High", "Low", "Close", "Volume", "Adjusted")

head(df)
#       Date   Open   High    Low  Close   Volume  Adjusted
# 1  2015/6/3 1322.0 1350.0 1319.5 1319.5 16815000   1319.5
# 2  2015/6/2 1313.0 1326.5 1300.5 1326.0 14139600   1326.0
# 3  2015/6/1 1293.0 1304.0 1282.5 1302.5  8270000   1302.5
# 4 2015/5/29 1300.0 1308.5 1289.5 1300.5 10061300   1300.5
# 5 2015/5/28 1298.5 1312.0 1293.0 1299.0 14069200   1299.0
# 6 2015/5/27 1270.0 1291.0 1262.0 1280.5 12494400   1280.5

candle_stick <- read.zoo(df, header = T)

par(family = "HiraKakuProN-W3")
candleChart(candle_stick, name = "日産自動車", multi.col = TRUE, theme = "white")

7201.T

時系列解析

時系列解析 (Time Series Analysis) は, ある現象の変動を過去の動きとの関連で捉えようとするものです。

まず, 時系列解析を行う上で, 時系列データが 定常過程 か 非定常過程 かがポイントとなります。

  • 定常 (Stationarity): 一見したところ不規則な現象でも時間的に変化しない一定の確率モデルの実現値とみなせるような平均や分散が時間変化しない時系列。白色雑音を含む。
  • 非定常 (Non-Stationarity): 平均の周りの変動の仕方が時間的に変化している時系列。平均や分散の時間変化がある場合。

つまり, 定常は平均、分散が時間変化せずに, 自己相関が時間差のみに依存する性質をさします。
また, 定常性は長期的には平均への回帰を期待しているとも言える。

tseries::adf.test() で単位根検定を行います。

df.zoo <- zoo::zoo(df$Close, order.by = as.Date(df$Date, format = '%Y/%m/%d'))
# to time-series
close <- ts(df.zoo)

head(close)

tseries::adf.test(close)
# Augmented Dickey-Fuller Test
#
# data:  close
# Dickey-Fuller = -2.6838, Lag order = 10, p-value = 0.2889
# alternative hypothesis: stationary

原系列は非定常時系列のようです。原系列が非定常の場合でも, そのd階差は定常となる場合があります。

close.diff <- diff(log(close))

tseries::adf.test(close.diff)
# Augmented Dickey-Fuller Test
#
# data:  close.diff
# Dickey-Fuller = -11.138, Lag order = 10, p-value = 0.01
# alternative hypothesis: stationary

p値が 0.01 であることから差分系列は定常時系列と考えられます。

時系列データのモデリング

今回のような株価や為替のデータの場合, あまりに古いデータから使ってしまうと経済の仕組みが変化していることがあり返って予測精度を落とす場合があります。
今回は 2010年以降のデータでモデリングしてみます。

今回試すモデルは以下です。

  • AR (Autoregressive): 自己回帰モデル
  • ARMA (Autoregressive and Moving Average): 自己回帰移動平均モデル
  • ARIMA (Autoregressive, Integrated and Moving Average): 自己回帰和分移動平均モデル

ARモデル

ARモデルは過去の観測値から定まるモデルで, 自己回帰の次数 m と自己回帰係数 a として下記で表されます。

     \begin{eqnarray*}   y_n = \sum_{i=1}^m a_i y_{n-1} + v_n \end{eqnarray*}

stats::ar() は特に指定しない限り ラグ次数 を AIC が最小となるように求めます。

ARMAモデル

ARモデルに移動平均 MA の項を追加したのがARMAモデルで, 移動平均の次数 l, 移動平均係数 b として下記となります。

     \begin{eqnarray*}   y_n = \sum_{i=1}^m a_i y_{n-1} + v_n - \sum_{i=1}^l b_i v_{n-1} \end{eqnarray*}

ARMAモデルの動機のひとつは AR で大きくなりすぎる次数を減らすことで, 移動平均を合成することで滑らかになります。

ARMAモデルは stats::arima() で作れます。

ARIMAモデル

ARIMAは AR (自己回帰), Integrated (和分), MA (移動平均) から成り, ARMAモデルに和分の I が追加されています。
和分とは, 定常な時系列を得るために時系列の観測値の階差をプロセスを指します。原系列 d-1 は非定常過程であっても, d階差分をとった過程が定常過程となる場合を d次和分過程 と言います。

ARIMA は近似的に定常な時系列が得られるまでd階差分をとり, その時系列に ARMA を当てはめます。つまり, d次和分過程 にARMAモデルを適用したものになります。ARIMA(m,0,l) の場合 ARMA(m,l) と一致します。

require(forecast)
# Plots a time series along with its acf (AutoCorrelation Function) and either its pacf (Partial ACF),
# lagged scatterplot or spectrum.
par(family="HiraKakuProN-W3")
tsdisplay(data.ts, main = "日産")

model <- auto.arima(
    data.ts,
    max.p = 2, # AR degree
    max.q = 2, # MA degree
    ic = "aic",# model selection.
    trace = T, # reported option
    stepwise = T
)

summary(model)

forecast::tsdisplay() でACF, PACFをプロットします。
ACF (自己相関関数)は, 自身のラグの時系列との相関, PACF (偏自己相関係数)は, 時系列とそれ以前のラグでは説明できない相関相関量です。

ACF

forecast::auto.arima()でARIMAモデルのモデル選択を自動で行います。icで基準を指定します。AIC基準の場合は下記です。

model <- auto.arima(
    df.ts,
    max.p = 2, # AR degree
    max.q = 2, # MA degree
    ic = "aic",# model selection.
    trace = T, # reported option
    stepwise = T
)

summary(model)
# ARIMA(2,1,2) with drift         : 11139.28
# ARIMA(0,1,0) with drift         : 11134.47
# ARIMA(1,1,0) with drift         : 11136.04
# ARIMA(0,1,1) with drift         : 11135.78
# ARIMA(0,1,0)                    : 11133.21
# ARIMA(1,1,1) with drift         : 11137.57
#
# Best model: ARIMA(0,1,0) with drift

ARIMA(0,1,0) が選択されました。これは, 1階差分の AR(0), MA(0) です。

sigma^2 estimated as 260.4:  log likelihood=-5136.17
AIC=10276.34   AICc=10276.35   BIC=10286.56

Training set error measures:
                       ME     RMSE      MAE         MPE     MAPE     MASE
Training set 0.0005219955 16.12392 12.06049 -0.02388146 1.412625 1.000159
                   ACF1
Training set 0.02431441
Training set 0.09790725

forecast() で予測を行います。levelで信頼区間の幅を指定, hで予測する期間を指定できます。

f <- forecast(
    model,
    level = c(50, 95), # confidence interval
    h = 20 * 3 # after 3 months
)
plot(f)

ARIMA

stl() で原系列を 季節変動 + トレンド + 残差 に分解して図示してみます。

# data = seasonal + trend + residualError
stl <- stl(data.ts, "periodic")
plot(stl)

trend

CodeはGitHubに置いてます。

今回試した AR, ARMA, ARIMA モデルは単変量を扱うモデルですが, この他に単変量を扱うモデルには GARCH があります。
一方, ARIMA は時系列の挙動あるいは変動要因に対して解釈を与えるのに難があるという見方もあり, モデリングの柔軟性と増減要因の説明力という点において状態空間モデルに利点があります。状態空間モデルは観測方程式と状態方程式を組み合わせたモデルで, 状態空間モデルはARIMAを内包します。

参考にさせて頂いた本は, 以下の3冊です。

時系列解析入門。

経済・ファイナンスデータの計量時系列分析。

ファイナンスのためのRプログラミング。


[1] Rで計量時系列分析:AR, MA, ARMA, ARIMAモデル, 予測
[2] Rで計量時系列分析:単位根過程、見せかけの回帰、共和分、ベクトル誤差修正モデル

1件のコメント

  1. ピンバック: R Data Mining の事例リンク | sleepingDog – please do not wake –

コメントを残す

必須欄は * がついています