【事例】SARIMAモデルで読書支出を予測してみた

RでのSARIMAモデルを使った実装例がWEB上にあまりなかったので投稿。
SARIMAモデル自体の説明は省きますが、心優しきかたが「[TV朝日の視聴率推移をSARIMAモデルで予測してみる](https://qiita.com/mshinoda88/items/749131478bfefc9bf365)」にSARIMAモデルの概要とPythonでの実行例をまとめてくださっているので、モデルの概要理解がまだの人はそちらもご参照ください。



# 対象読者
SARIMAモデル（季節性自己回帰移動平均モデル）で何か分析したいが、Rでの実行方法がわからない人。
SARIMAモデルで実際どの程度予測できるのかに興味がある人。


# データ準備
元データ：総務省統計局の家計調査（家計収支編）
時系列データ（二人以上の世帯）品目別の支出金額（2001年1月〜）を使用しました。
  http://www.stat.go.jp/data/kakei/longtime/csv/h-mon-a.csv
  
上記の元データは形式としてはcsvだが、かなりエクセルに最適化されて作ってあるため、そのままではRで処理するのが大変。なので、Rで読み込めるようにエクセル上で加工したデータを用います。

加工データ：下記、github上に配置しています。
https://github.com/chaplin96/household_expenses

加工データをダウンロードの上、下記コマンドを実行。

```R:script.R
# 今回使用するライブラリをインストールしておく。
install.packages("tseries")
install.packages("forecast")
install.packages("ggplot2")
install.packages("ggfortify")

library(tseries)
library(forecast)
library(ggplot2)
library(ggfortify)

# 家計支出データをデータフレーム型として読み込む
df_expenses<-read.csv("household_expenses.csv",fileEncoding = "Shift-JIS",stringsAsFactors = FALSE)

# tsオブジェクトに変換
ts_expenses<-ts(select(df_expenses,-time),start=c(2000,1),freq=12)

```


# データ観察


```R:script.R
# どのような品目があるかを観察
str(df_expenses)

# 実行結果
'data.frame':	217 obs. of  186 variables:
 $ time                                : chr  "2000-01-01" "2000-02-01" "2000-03-01" "2000-04-01" ...
 $ 世帯数分布.抽出率調整.              : int  10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 ...
 $ 集計世帯数                          : int  7887 7942 7934 7922 7928 7917 7907 7908 7917 7937 ...
 $ 世帯人員.人.                        : num  3.32 3.32 3.32 3.32 3.31 3.31 3.31 3.31 3.31 3.3 ...
 $ X18歳未満人員.人.                   : num  0.74 0.75 0.75 0.75 0.75 0.74 0.74 0.74 0.74 0.73 ...
 $ X65歳以上人員.人.                   : num  0.52 0.53 0.53 0.52 0.52 0.53 0.54 0.53 0.53 0.54 ...
 $ うち無職者人員.人.                  : num  0.41 0.41 0.41 0.41 0.41 0.42 0.42 0.42 0.42 0.42 ...
 $ 有業人員.人.                        : num  1.51 1.51 1.51 1.52 1.53 1.53 1.52 1.52 1.51 1.52 ...
 $ 世帯主の年齢.歳.                    : num  52.4 52.6 52.7 52.6 52.7 52.6 52.7 52.9 52.8 53 ...
 $ 持家率...                           : num  76 76.3 76.2 75.8 76.1 76 75.7 75.9 75.8 76.6 ...
 $ 家賃.地代を支払っている世帯の割合...: num  22.1 22 22.6 22.2 22.6 22.9 23 22.4 22.6 22.8 ...
 $ 消費支出                            : int  309621 290663 335341 335276 308566 297648 326480 309993 296457 309193 ...
 $ 食料                                : int  73580 73309 79726 77344 81415 76721 83782 86988 76985 79477 ...
 $ 穀類                                : int  6100 6915 7496 7470 7447 7328 7794 7324 7631 8314 ...
 $ 米                                  : int  2338 2919 3226 3346 3264 3165 3160 3164 3867 4386 ...
・・・(以下省略）


# 例として、今回は読書支出を取り上げる。
books<-ts_expenses[,"読書"]

autoplot(
  books,
  main="読書支出",
  xlab="year",
  ylab="expenses"
)

```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254014/ce3f5956-070f-9cd8-1006-ff9bcd0dcdcf.png)


```R:script.R
# 月別でグラフを描く
ggsubseriesplot(books)

# 一階差分をとったデータのグラフを描く
diff_books<-diff(books,lag=1)
ggtsdisplay(diff_books)


# 季節階差をとったデータのグラフを描く
seas_diff_books<-diff(books,lag=frequency(books))
ggtsdisplay(seas_diff_books)

```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254014/c1a73795-4a2b-f53b-e35b-a46798bd22ad.png)

月別グラフから読み取れること
(1)1月に最も少なく、12月に最も多くなる
(2)どの月においても、年々減少傾向である
  
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254014/83d83071-fa2c-4e83-6390-94f59b0f522a.png)

一階差分をとったデータのグラフから読み取れること
(1)一階差分をとると定常過程になっている
(2)数期前のデータと自己相関がある
(3)12期前（1年前）のデータと強い自己相関がある


![image.png](https://qiita-image-store.s3.amazonaws.com/0/254014/803e4a18-24fe-0611-dfc2-903b74b1bdc6.png)

季節階差をとったデータのグラフから読み取れること
(1)1~5期前のデータと自己相関がある
(2)季節階差をとることで減少したが、相変わらず１２期前（1年前）のデータと自己相関がある
(3)原系列の減少トレンドを、季節階差の定数項として表すことができそう

# SARIMAモデル構築

```R:script.R
# 学習データとテストデータに分ける
train_books<-window(books,end=c(2015,12))
test_books<-window(books,start=c(2016,1),end=c(2017,12))


# forecastパッケージに実装されているauto.arima関数を用いてSARIMAモデルを構築する。AIC(又はBIC)が最小となる次数のモデルを採用する。
m_sarima_books<-auto.arima(
  y = train_books,
  ic = "aic",
  D = 1,
  max.order = 5,
  seasonal = T,
  stepwise=F,
)
# auto.arima関数の引数の説明
# ic:どの情報量基準を使うかを指定(aicの他に、bicがある)
# D：何度季節階差をとるかを指定する（原系列とseas_diff_booksを見比べて、1階季節階差をとった方が、減少トレンドをうまく表現できそうだっため、今回は1を指定してみた）
# max.order:ARモデルの次数p、MAモデルの次数q、SARモデルの次数P、SMAモデルの次数Qの和の最大値を指定
# seasnal:季節性を考慮するかを指定。FALSEを指定すると、季節成分を考慮しない（P,D,Qを用いない)モデルを構築する
# stepwise:次数選択において、ステップワイズ法を使うか。falseならステップワイズを用いずに、全ての組み合わせに関して計算する

m_sarima_books
# 実行結果
Series: train_books 
ARIMA(0,1,1)(2,1,1)[12] 

Coefficients:
          ma1     sar1     sar2     sma1
      -0.6791  -0.1466  -0.1796  -0.6052
s.e.   0.0609   0.1274   0.1020   0.1210

sigma^2 estimated as 7702:  log likelihood=-1058.17
AIC=2126.35   AICc=2126.69   BIC=2142.28


# 残差に自己相関が残っていないか、正規分布しているかをチェック
checkresiduals(m_sarima_books)

```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254014/c2d7ef76-40c8-70d0-60e0-5801a813d3fe.png)

６期ほど前に自己相関が少し残っているのが気になりますが、今回は許容範囲として先に進みます。


# 予測 
```R:script.R
# 構築したモデルを使って、2年(24期)先まで予測してみる
f_sarima_books<-forecast(
  m_sarima_books,
  h=24,
  level = c(95,70)
)
# levelの95,75は95%,75%予測区間も併せて出力するという指定

# 予測と実際の値をグラフで比較してみる。青の点線が予測値。
plot(window(books,end=c(2017,12)),xlab="時間",ylab="読書支出(円)",main="予測と実際の比較")
lines(f_sarima_books$mean,col="blue",lwd=2,lty="dotted")

```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254014/6640f1bc-f12e-540c-8874-597a94f1cfbc.png)

ある程度うまく予測できてそうですが、次の節ではどれくらいうまく予測できたかを数値で表現します。

# 精度検証

```R:script.R
# forecastパッケージのaccuracy関数を用いて予測精度を評価する
accuracy(f_sarima_books,test_books)

# 実行結果
                    ME     RMSE      MAE         MPE     MAPE      MASE         ACF1 Theil's U
Training set  -3.663173 83.78706 60.08917 -0.09521795 1.394424 0.5650724 -0.005997121        NA
Test set     -43.082400 78.31399 61.46071 -1.22261005 1.741044 0.5779702 -0.044425539 0.4378374
```

様々な指標が出てくるが、元データと単位が揃っているRMSE（平均二乗誤差平方根）が直感的にはわかりやすいと思います。

月ごとの読書支出を80円程度の誤差で予測できているので、精度は良さそうです。


しかし、今回構築したSARIMAモデルが本当に優れているのかは、他の単純なモデルと比較して初めてわかります。
例えば、今までの月平均を予測値として出力するモデルのRMSEと、SARIMAのRMSE比べてみることにしましょう。

```R:script.R
# 月毎の平均値を算出する
f_mean<-tapply(train_books,cycle(train_books),mean)

# 2年分の予測値を時系列データオブジェクトとして作成する
f_mean_ts<-ts(rep(f_mean,2),start=c(2016,1),frequency = 12)

# 予測精度の検証
accuracy(f_mean_ts,test_books)
# 実行結果
               ME     RMSE     MAE       MPE     MAPE      ACF1 Theil's U
Test set -732.875 738.9992 732.875 -20.71808 20.71808 0.3346646  4.019559

```
単純な平均値では、RMSEは740程度なので、構築したSARIMAモデルは比較的良い精度だと言えそうです。（めでたし。）

# まとめ
SARIMAモデルの構築自体は、forecastパッケージを用いれば簡単にできます。ただし、精度の良いモデルを作るには色々な角度から作図してデータの特性を理解する必要があります。

また、今回使用したデータの中には読書支出以外にも多くの支出データが含まれているので、ぜひ遊んでみてください！！
余談ですが、VAR（ベクトル自己回帰）でデータ同士の因果関係を分析してみても面白いかもしれませんね。


# 参考文献
馬場真哉『時系列分析と状態空間モデルの基礎 RとStanで学ぶ理論と実装』
沖本竜義『経済・ファイナンスデータの軽量時系列分析』
「[TV朝日の視聴率推移をSARIMAモデルで予測してみる]
(https://qiita.com/mshinoda88/items/749131478bfefc9bf365)」