Rでar関数や残差分析などの復習

R Python 時系列分析 統計 統計分析

前回「Rのts.arima関数を使ってみる」を書いてから、5日経過。
すでに色んなことを忘れているｗｗｗ。1日開けば殆ど忘れます。SE的にすべきことが多いので。

今回は、まずRの使い方の復習から。

>setwd("D:\\dev\\DataAnalysis\\Data")
>f1min_20151112_Day<-read.csv("f1min_20151112_Day.csv")
>f1min_20151112_Day.close<-f1min_20151112_Day[,5]
>f1min_20151112_Day.close.ts<-as.ts(as.numeric(f1min_20151112_Day.close))
>library(forecast)
>f1min_20151112_Day.close.ts.arima<-auto.arima(f1min_20151112_Day.close.ts,max.p=100,stepwise=T,trace=T)
>plot(forecast(auto.arima(f1min_20151112_Day.close.ts,max.p=100,stepwise=T,trace=T),level=c(10,50,90),h=100))

とりあえず、Rでグラフ描画。ここまでは前回と同じ（ファイル名など若干変更しているがデータは同じ）。

◎自己相関関数（ACF：Auto Correlation Function）の描画

> acf(f1min_20151112_Day.close.ts)

離れたラグでも自己相関がある？

◎偏自己相関関数(PACF：partial autocorrelation function)の描画

> pacf(f1min_20151112_Day.close.ts)

直前のラグ以外は相関が非常に小さいので、直前データからの影響のみ考えれば良い？

◎モデル推定
試しにARモデルを仮定する。

Rの関数ARで自己回帰モデルを求めてみる。

> ar(f1min_20151112_Day.close.ts,aic=TRUE, order.max=25)

Call:
ar(x = f1min_20151112_Day.close.ts, aic = TRUE, order.max = 25)

Coefficients:
      1        2        3        4        5        6        7  
 0.7873   0.1407   0.1064  -0.0722  -0.0785  -0.0521   0.1434  

Order selected 7  sigma^2 estimated as  65.73

次数7のARモデルが提案された。

AICを確認してみる。

> ar(f1min_20151112_Day.close.ts,aic=TRUE, order.max=25)$aic
          0           1           2           3           4           5 
896.2693150   8.2582626   0.6746665   1.5664388   2.4096834   4.3744895 
          6           7           8           9          10          11 
  5.0628233   0.0000000   1.9851290   3.4629374   3.7822072   5.7391055 
         12          13          14          15          16          17 
  5.4711945   5.8728777   7.1771941   7.5596980   9.5242780  10.5502276 
         18          19          20          21          22          23 
 10.5247654  11.2311400  12.4292297  14.0319148  15.9847824  17.6573267 
         24          25 
 17.6592898  19.2360679

プロットしてみる。

> plot(ar(f1min_20151112_Day.close.ts,aic=TRUE, order.max=25)$aic)


次数0だけが非常にAICが大きいので、これを除いてプロットしてみる。

> plot(ar(f1min_20151112_Day.close.ts,aic=TRUE, order.max=25)$aic[2:24])


次数2のARモデルが次数7のARモデルの次にAICが小さいようだ。

◎モデルの診断
残差分析を行って、作成したモデルの適切さを判断する。
ARモデルは、「残差が平均0の正規分布に従い、互いに独立であること」を条件としたモデルなので、これに近いことが望まれる。

リュング・ボックス検定を行って、残差の独立性を検定してみる。

> f1min_20151112_Day.close.ar<- ar(f1min_20151112_Day.close.ts,aic=TRUE,order.max=25)
> Box.test(f1min_20151112_Day.close.ar$res, type="Ljung")

        Box-Ljung test

data:  f1min_20151112_Day.close.ar$res
X-squared = 0.094802, df = 1, p-value = 0.7582

ここで以下のwikipediaからの引用を参照する。
これを見ると、P値が大きく帰無仮説が棄却できないと、「残差は無作為である」ことが言えそう。
大きいとは定量的に幾ら以上のことを言うのかは。。。わからない。。。
有意水準をαとして、p値がα以上の場合に帰無仮説を棄却できない（逆に、p<αの際に帰無仮説が棄却される）のだが、ここでは有意水準αを0.05(5%)として、p値は大きいとしておこう・・・。

Ljung-Box 検定は、以下のとおり定義される。
H0：（帰無仮説）データは無作為である。
Ha：（対立仮説）テータは無作為でない。
ここで、Ljung-Box 統計量（-とうけいりょう、英: Ljung-Box statistic）を以下のとおり定義する。

これで、次数7のARモデルは残差の独立性についてはOKであろう。
残差の正規性を検証する。


さて、ar関数を使用してAR(7)を求めたが、arima関数を使用して同じものを求めてみる。
何故かと言うと、AICの計算値はAR(7)の次にAR(2)が大きかったので、AR(2)も調べてみたいのだが、ar関数を使用したAR(2)の求め方が分からないからｗ。なので、arima関数を使用してAR(7)とAR(2)を求めてみようということ。

> f1min_20151112_Day.close.arima700<- arima(f1min_20151112_Day.close.ts,orde=c(7,0,0))
> f1min_20151112_Day.close.arima700

Call:
arima(x = f1min_20151112_Day.close.ts, order = c(7, 0, 0))

Coefficients:
         ar1     ar2     ar3      ar4      ar5      ar6     ar7   intercept
      0.7778  0.1684  0.1130  -0.0840  -0.0589  -0.0934  0.1544  19684.1574
s.e.  0.0547  0.0696  0.0701   0.0703   0.0703   0.0700  0.0554     15.7499

sigma^2 estimated as 57.13:  log likelihood = -1171.65,  aic = 2361.3

これがAR(7)。ar関数を使用して求めた結果と若干数値が異なるが、気にしないｗ。

これに残差分析を行う。

> Box.test(f1min_20151112_Day.close.arima700$resid, type="Ljung")

        Box-Ljung test

data:  f1min_20151112_Day.close.arima700$resid
X-squared = 0.002245, df = 1, p-value = 0.9622

有意水準α=0.05として、p値は0.9622>αなので、帰無仮説を棄却できない。
やはり、arima関数を使用した場合でも残差の独立性に問題はない。

ちなみに、arima関数を使用した場合、モデル診断を行うツールとしてtsdiag関数を使用できる。

> tsdiag(f1min_20151112_Day.close.arima700, gof.lag=12)

◎予測
predict関数を使用すればできる。
n.aheadには予測期間を指定する。

> (f1min_20151112_Day.close.ar.predict<-predict(f1min_20151112_Day.close.ar,n.ahead=10))
$pred
Time Series:
Start = 341 
End = 350 
Frequency = 1 
 [1] 19686.72 19687.08 19685.90 19688.42 19690.72 19693.76 19692.16 19692.42
 [9] 19692.48 19691.69

$se
Time Series:
Start = 341 
End = 350 
Frequency = 1 
 [1]  8.107262 10.318167 12.019842 13.719601 15.042394 15.963852 16.562025
 [8] 17.229448 17.821954 18.355708

プロットしてみる。

> ts.plot(f1min_20151112_Day.close.ar.predict$pred)

過去の実測値に予測データを付け加えてプロットしてみる。

> ts.plot(f1min_20151112_Day.close)
> lines(f1min_20151112_Day.close.ar.predict$pred, col="red")