ARIMAXモデルにおいて、X(外部要因)に何を入れるか考える際にいちいち1つ1つ試していくのは手間がかかる。そのため、loop処理で一気に計算をおこない結果を比較したい。

データ

テキトーに以下のデータで考える。

> df
# A tibble: 7 x 5
  date       value    x1  x2_1  x2_2
  <date>     <dbl> <dbl> <dbl> <dbl>
1 2017-04-01  126.  2.82  850.   412
2 2017-05-01  132.  2.04  511.   249
3 2017-06-01  130.  2.28  698.   340
4 2017-07-01  122.  1.24  323.   157
5 2017-08-01  143.  4.56 1175.   574
6 2017-09-01  172.  1.69  425.   207
7 2017-10-01  140.  1.41  285.   134

事前準備

# xtsに変換
xts_df <- as.xts(read.zoo(df))

# testとtrainに分割
xts_df.train = xts_df[1:5]
xts_df.test = xts_df[6:7]

外部要因部分

ARIMAXモデルをauto.arimaで作成する際には、以下のようにxregオプションに外部要因となるデータを渡す必要がある。

# xregに渡したいデータ
x_data = xts_df.train[,c('x1', "x2_1")]

  #=>
  # x1      x2_1
  # 2017-04-01 2.821320  849.7471
  # 2017-05-01 2.038116  510.6461
  # 2017-06-01 2.283226  698.1490
  # 2017-07-01 1.235125  323.1432
  # 2017-08-01 4.555123 1174.8850
  # 2017-09-01 1.694707  424.8735

# ARIMAXモデルの自動作成
  model_arimax = auto.arima(
    y = xts_df.train[,'value'],
    xreg = x_data,
    ic = "aic",
    max.order = 20,
    stepwise = F,
    approximation = F,
    parallel = T,
    num.cores = 2
  )

今回は、このxregに渡したいデータ(x_data)をloop処理で変えていく。
変えていくx_dataは、欲しいカラム部分のベクトルが書かれたlistを作成し、そのlistをloopで変えていく。

# xregに渡すデータに入れたいカラム組み合わせを複数記載
models_list = list(model1 = c('x1'),
                   model2 = c('x1', "x2_1"),
                   model3 = c('x1', "x2_2")
)

for (v in models_list) {
  x_data = xts_df.train[,v]
  
  # ARIMAXモデルの自動作成
  model_arimax = auto.arima(
    y = xts_df.train[,'value'],
    xreg = x_data,
    ic = "aic",
    max.order = 20,
    stepwise = F,
    approximation = F,
    parallel = T,
    num.cores = 2
  )
}

これでloop毎に、model1, model2, model3で指定したカラムのデータが外部要因となったARIMAXモデルがそれぞれで作成される。

ただし、このままではそれぞれのloopでの結果は保存されないので結果を保存する必要がある。

結果の保存

構築したモデルは上記コードのようにmodel_arimaxオブジェクトに格納されているとする。その場合、auto.arimaの結果は以下のようにして取り出すことができる。

  # 結果
  arma =  model_arimax$arma #次数の結果(前から順に AR次数，MA次数，季節要素のAR次数，季節要素のMA次数，季節要素の期間，差分，季節要素の差分)
  aic = model_arimax$aic # AIC
  coef = model_arimax$coef #係数推定値

[次数の結果の参考]

auto.arimaで学習した結果はあるんだけど，best modelの次数を忘れちゃった時 - Qiita

また、構築されたモデルをtestデータで予測した値に対するRMSEは以下のようにして取り出すことができる。

  # testで予測
  model_arimax_test = forecast(
    model_arimax,
    xreg = xts_df.test[,v],
    h = 1,
    level = c(95,70)
  )

# testで予測した値に対するRMSE
rmse = accuracy(model_arimax_test, xts_df.test[,'value'])['Test set','RMSE']

ちなみにRMSEの意味は以下参考。

いろいろな誤差の意味（RMSE、MAEなど） - 具体例で学ぶ数学

結果の保存(いい感じの保存）

それぞれの結果は上記コードから取り出せるが、見づらいのでtibble形式で見やすいようにデータを加工して保存する。

  # モデルの推定値をtidyに整形
  temp_model_coef = tidy(model_arimax) 
  
  # 結果を整形して保存
  model = list(model = temp_model_coef) %>%
    enframe(name = 'model_no', value = 'coef') %>% 
    mutate(model = paste(v, collapse = ' + '), #モデル外部要因
           # testで予測した結果に対するRMSE
           rmse = accuracy(model_arimax_test, xts_merged.test[,'recall_townwork'])['Test set','RMSE']
    )

全容

今までの内容をまとめると以下のようになる。

xts_df <- as.xts(read.zoo(df))

# testとtrainに分割
xts_df.train = xts_df[1:5]
xts_df.test = xts_df[6:7]

# 入れたい外部要因
models_list = list(model1 = c('x1'),
                   model2 = c('x1', "x2_1"),
                   model3 = c('x1', "x2_2")
)

n = 0
models = tibble()

for (v in models_list) {
  print(v)
  n = n + 1
  
  x_data = xts_df.train[,v]
  
  
  # ARIMAX
  model_arimax = auto.arima(
    y = xts_df.train[,'value'],
    xreg = x_data,
    ic = "aic",
    max.order = 20,
    stepwise = F,
    approximation = F,
    parallel = T,
    num.cores = 2
  )
  
  # 結果
  arma =  model_arimax$arma #次数の結果(前から順に AR次数，MA次数，季節要素のAR次数，季節要素のMA次数，季節要素の期間，差分，季節要素の差分)
  aic = model_arimax$aic # AIC
  coef = model_arimax$coef #係数推定値 
  
  # testで予測
  model_arimax_test = forecast(
    model_arimax,
    xreg = xts_df.test[,v],
    h = 1,
    level = c(95,70)
  )
  
  # tidyに整形
  temp_model_coef = tidy(model_arimax) 
  
  # 欲しいデータを整形して保存
  model = list(model = temp_model_coef) %>%
    enframe(name = 'model_no', value = 'coef') %>% 
    mutate(model_no = n,
           model = paste(v, collapse = ' + '), #モデル外部要因
           # testで予測した結果に対するRMSE
           rmse = accuracy(model_arimax_test, xts_df.test[,'value'])['Test set','RMSE'],
         aic = aic,
         arma = paste(arma,collapse = ',')
    )
  
  # 各loop結果を保存
  models = rbind(models, model)
}

結果の見方

modelsオブジェクトに各結果が保存された。なお、coef列は各モデルの係数が別途入っている

> models
# A tibble: 3 x 6
# model_no coef             model      rmse   aic arma         
# <dbl> <list>           <chr>     <dbl> <dbl> <chr>        
#   1        1 <tibble [3 × 3]> x1         34.7  32.3 1,0,0,0,1,0,0
# 2        2 <tibble [3 × 3]> x1 + x2_1  32.2  31.5 0,0,0,0,1,0,0
# 3        3 <tibble [3 × 3]> x1 + x2_2  32.1  31.7 0,0,0,0,1,0,0

coefの中身を見たい場合はunnestすると、nestされたcoefが展開された状態で表示される。

df_result = models %>% 
  dplyr::select(model_no, coef) %>% 
  unnest %>% 
  mutate_if(is.double, round, digits=2)  #下二桁表示に変換


> df_result
# A tibble: 9 x 4
# model_no term      estimate std.error
# <dbl> <chr>        <dbl>     <dbl>
#   1        2 ar1          -0.77      0.3 
# 2        2 intercept   117.        3.33
# 3        2 x1            5.19      1.38
# 4        3 intercept   118.        3.01
# 5        3 x1           15.8       5.07
# 6        3 x2_1         -0.04      0.02
# 7        4 intercept   117.        3.05
# 8        4 x1           15.9       5.31
# 9        4 x2_2         -0.08      0.04

ただし、上記では見づらいので以下のようにspreadで各係数に対して、モデル毎に横持ちにすると見やすくなる。

df_result = models %>% 
  dplyr::select(model_no, coef) %>% 
  unnest %>% 
  mutate_if(is.double, round, digits=2) 

df_result %>% 
  mutate(term = fct_inorder(term)) %>%  #デフォルトではアルファベット順なので出てきた順にする
  dplyr::select(model_no, term, estimate) %>% 
  spread(model_no, estimate)

# A tibble: 5 x 4
# term          `2`     `3`     `4`
# <fct>       <dbl>   <dbl>   <dbl>
#   1 ar1         -0.77   NA      NA   
# 2 intercept  117.    118.    117.  
# 3 x1           5.19   15.8    15.9 
# 4 x2_1        NA      -0.04   NA   
# 5 x2_2        NA      NA      -0.08