モモノキ＆ナノネと一緒にRで時系列分析を試してみよう
Rで時系列分析、AirPassengersデータとSARIMAモデルを使って予測をやってみよう
ナノネ、統計解析用フリーソフト『R』の使い方を練習するよ。今回は時系列分析をやってみよう。
時系列分析だね！データはPythonでも使ったAirPassenger？
うん。Rだとデフォルトでデータが入ってるから読み込みもラクチンだね。
さっそくAirPassengersのデータをチェックしてみよう。
AirPassengers
飛行機乗客数のデータで1949年から1960年までの月次データ。数値は人数で、たしか千人単位だったよね。
plot(AirPassengers)
Rにおまかせでグラフにプロットしてみた。
class(AirPassengers)
時系列データは『ts』とうい形式になるみたいだよ。
str(AirPassengers)
tsはTime-Seriesの略みたいだね。
Rで時系列分析をやるのは今回が初めてだから、まだやりかたが良くわからないけど．．．。グーグル先生に質問しながら、とりあえず実践してみよう。
ラジャー。
まずはデータ準備。AirPassengersは144カ月、12年分のデータがあるね。モデルで予測もやりたいから、モデル作成用のデータは最後の1年分を除外しておこう。
data <- AirPassengers[1:132] # 1:(144-12)
最後の1年分を除外して、こんな感じかな？？
print(data)
class(data)
モモノキ、ダメみたい。ts形式だったのに、ただのベクトルデータになっちゃった。
時系列データはts()を使ってデータを加工するといいみたい。
data <- ts(AirPassengers[1:132])
data
OK、ts形式に戻った。でも年月の設定が消えちゃったけど？
ts()のオプションで設定できるよ。startとfrequencyに値を渡せばいいみたい。今回の場合はstartにベクトルで開始の年と月、frequencyは月次データだから12を設定すればOK。
data <- ts(AirPassengers[1:132], start=c(1949, 1),frequency=12)
data
うまくいったかも。
次はモデルの作成だけど、そのまえに準備。『forecast』というのを使うからlibrary()で読み込んでおいてね。
# install.packages("forecast") # 必要に応じてインストール
library(forecast)
forecast準備OK。
参考までにstl()を使うと季節変動、トレンド、残差にデータを分離できるよ。グラフにも簡単にプロットできて、上から順に生データ、季節変動、トレンド、残差だよ。データチェックに活用できるね。
data.stl <- stl(data,s.window="periodic")
plot(data.stl)
Pythonで出力したグラフと一緒だ。
さっそくモデルをつくるよ。Rでも(S)ARIMAモデルを使ってみよう。
今回もまたモデルの最適次数とか探すの？Pythonで試したときけっこう手間取ったよね。
auto.arima()を実行すれば推奨次数をRが全部を教えてくれるから簡単だよ。auto.arimaのパラメータはコメントを参考にしてね。
model.arima <- auto.arima(
data, # データ指定
ic = "aic", # 情報量基準の指定
trace = T, # トレース出力（結果一覧）
stepwise=T, # 計算時の次数パターン省略
approximation=F, # 近似計算の適用
seasonal=T # 季節調整
)
ARIMA(1,1,0)(0,1,0)[12]がベストだって言ってる。名前はARIMAだけど季節調整(SARIMA）にも対応してるんだ。
うん、季節調整の次数チェックもしてくれるから便利だね。もし手動でパレメータを設定したい場合は、下みたいなコードでモデルをつくれるよ。
arima(
data,
order=c(1,1,0),
seasonal=list(order=c(0,1,0))
)
tsdiag()で作成したモデルの残差などをグラフ化できるよ。
tsdiag(model.arima)
グラフが3個あるけど、見かたは？
一番上が残差プロット、2番目が残差の自己相関、3番目がリュング・ボックス検定だよ。
残差と残差の自己相関は小さいほどいいよ。検定はホワイトノイズか？のチェック的なもので、白点が青破線を下回らなければとりあえずOKだよ。
次は作成した(S)ARIMAモデルで予測をやってみよう。forecast()で予測を実行できるよ。引数はモデルとlevelに信頼区間、hに予測期数を指定すればOK。
pred = forecast(model.arima, level=95, h=12)
pred
class(pred)
そのままグラフにプロットもできるよ。
plot(pred)
Rだと信頼区間も簡単にプロットできて便利だね。
後半期間をグラフで拡大、予測期間もあと1年増やして2年分でやってみよう。
plot(
forecast(model.arima, level = c(95), h = 24),
flty = 2,
xlim = c(1958, 1962),
xaxp = c(1958, 1962, 4),
ylim = c(200, 800),
yaxp = c(200, 800, 3),
)
lines(AirPassengers, lwd=2)
legend("topleft",
legend = c("pred", "actual"),
lty = c(2, 1),
lwd = 2,
col = c(4, 1),
cex = 1.5,
bty = 'n',
inset = c(0.05, 0.03)
)
けっこういい感じ？に予測できたかも。
今回のRで時系列分析は以上で。
またね！
0 件のコメント :
コメントを投稿