【統計】時系列分析（パネルデータ分析）その②

R、R Studio 統計

パネルデータ分析のメモ。その２

【目次】


以下で基本部分を取りまとめた、今回は実データで推定してみる。
【統計】時系列分析（パネルデータ分析） - こちにぃるの日記

使用するデータ

Rのplmパッケージの Grunfeld を用いる。
目的変数はinv、説明はvalue captal、属性情報はfirm、経時情報はyearに含まれる。
・firm　：企業（ファーム）を表す変数。
・year　：年度。
・inv　　：企業ごとの投資額。
・value　：企業の市場価値を表す変数。
・capital：企業の設備や設備の在庫を表す変数。

10社の1935年～1954年のデータが含まれる。
1950年までのデータを使ってパネルデータ分析をしてみる。

library(plm)
data("Grunfeld", package="plm")

Grunfeld$FLG <- ifelse(Grunfeld$year<=1950,1,0)
ADTS_T <- pdata.frame(Grunfeld[Grunfeld$FLG==1,], index=c("firm","year"))


パネルデータ分析の実装

固定効果モデルの実行。
valueとcapitalの両方が統計的に有意である。
価値valueが1単位増加すると、投資 inv は約0.0752単位増加する
資本capitalが1単位増加すると、投資 inv は約0.2063単位増加する

> #-- 固定効果モデル
> fit1 <- plm(inv ~ value + capital, data=ADTS_T, model="within")
> summary(fit1)
Oneway (individual) effect Within Model

Call:
plm(formula = inv ~ value + capital, data = ADTS_T, model = "within")

Balanced Panel: n = 10, T = 16, N = 160

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-154.17714  -11.85856    0.15505   10.77434  142.68974 

Coefficients:
        Estimate Std. Error t-value  Pr(>|t|)    
value   0.073680   0.010987  6.7061 3.945e-10 ***
capital 0.204658   0.022589  9.0599 7.158e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    377500
Residual Sum of Squares: 222510
R-Squared:      0.41058
Adj. R-Squared: 0.36677
F-statistic: 51.5477 on 2 and 148 DF, p-value: < 2.22e-16


変量効果モデルの実行。
こちらも同様にvalueとcapitalの両方が統計的に有意である。
価値valueが1単位増加すると、投資 inv は約0.0813単位増加する
資本capitalが1単位増加すると、投資 inv は約0.2083単位増加する

> #-- 変量効果モデル
> fit2 <- plm(inv ~ value + capital, data=ADTS_T, model="random")
> summary(fit2)
Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = inv ~ value + capital, data = ADTS_T, model = "random")

Balanced Panel: n = 10, T = 16, N = 160

Effects:
                  var std.dev share
idiosyncratic 1503.42   38.77 0.212
individual    5596.92   74.81 0.788
theta: 0.8715

Residuals:
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-133.3387  -13.4845   -0.5492    8.9296  158.9204 

Coefficients:
              Estimate Std. Error z-value Pr(>|z|)    
(Intercept) -5.1914084 26.5668704 -0.1954   0.8451    
value        0.0813006  0.0095707  8.4948   <2e-16 ***
capital      0.2082550  0.0221473  9.4032   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    441150
Residual Sum of Squares: 237130
R-Squared:      0.46247
Adj. R-Squared: 0.45562
Chisq: 135.077 on 2 DF, p-value: < 2.22e-16


ハウスマン検定の実行。
今回の場合、p値=0.1931で、帰無仮説が棄却されないため、変量効果モデルが妥当であると判断される。
ただし、ハウスマン検定や変量効果モデルは注意が多い。現実的に相関がない仮定は考えづらく、固定効果モデルを使う方が一般的？

> #-- ハウスマン検定
> phtest(fit1, fit2)

    Hausman Test

data:  inv ~ value + capital
chisq = 3.289, df = 2, p-value = 0.1931
alternative hypothesis: one model is inconsistent


推定モデルでの予測

1950年までのデータで作成したモデルで、1951年以降を予測してみる。
事前準備として、pdata.frameを使って専用のデータフレームを作成する。index に属性情報と経時情報を指定する。
作成したデータフレームをもとにpredictを使って予測値を算出する。

> #-- 評価指標
> ADTS_P <- pdata.frame(Grunfeld[Grunfeld$FLG==0,], index=c("firm","year"))
> pred1 <- predict(fit1, newdata=ADTS_P)
> pred2 <- predict(fit2, newdata=ADTS_P)


予測値と実測値でのRMSE、MAE、MAPEを算出。
外れ値の考慮や実測値に0を含むかでどの指標を使うか検討。

> ADTS_P$pred1 <- pred1
> ADTS_P$pred2 <- pred2
> 
> fnc1 <- function(OBS,PRED){
+   RMSE <- sqrt(1/length(PRED) * sum(OBS - PRED)^2)
+   MAE  <- 1/length(PRED) * sum(abs(OBS - PRED))
+   MAPE <- 1/length(PRED) * sum(abs(OBS - PRED) / PRED) * 100
+   print(paste0("RMSE:",round(RMSE,4)," MAE:",round(MAE,4)," MAPE:",round(MAPE,4)))
+ }
> fnc1(ADTS_P$inv, ADTS_P$pred1)
[1] "RMSE:204.1619 MAE:61.3335 MAPE:27.1647"
> fnc1(ADTS_P$inv, ADTS_P$pred2)
[1] "RMSE:181.0856 MAE:109.2446 MAPE:52.0501"


参考

統計学（出版：東京図書）, 日本統計学会編
https://www.jil.go.jp/institute/zassi/backnumber/2015/04/pdf/006-009.pdf

【統計】時系列分析（パネルデータ分析）

R、R Studio SAS 統計

パネルデータ分析のメモ。

【目次】

クロスセクション、時系列、パネルデータの違い

クロスセクションデータ（Cross Section Data）
　特定時点・複数変数のデータ。時間軸が無い。

時系列データ（Time Series Data）
　単一変数の経時的なデータ。

パネルデータ（Panel Data）
　複数変数の経時的なデータ。クロスセクションデータ＆時系列データ。



固定効果モデルと変量効果モデル

パネルデータに関する回帰モデルを次の通り表す。

　 y_{i,t}=\alpha + \beta X_{i,t} + u_{i,t}
　 u_{i,t}=\mu_{i} + \gamma_{i,t}

　 i：個体， t：時期， u_{i,t}：誤差項
　 X_{i,t}：説明変数， \mu_{i}：個別効果（Individual Effect）， \gamma_{i,t}：ランダム誤差

誤差項は「個別効果 \mu_{i} ＋ランダム誤差 \gamma_{i,t} 」で構成されると考える。
ここで、個別効果 \mu_{i} の相関の考えにて、固定効果モデルと変量効果モデルの２つのモデルがある。

固定効果モデル（Fixed Effect Model）
　個別効果 \mu_{i} と説明変数 X_{i,t} が相関する（内生性（Endogeneity）を持つ）と仮定する。
　個別効果 \mu_{i} は未知の定数。
　　⇒「未知の定数」の仮定のため、階差や個体ごとの時間平均を使って除去できる。

変量効果モデル（Random Effect Model）
　個別効果 \mu_{i} と説明変数 X_{i,t} は相関しない（独立である）と仮定する。
　個別効果 \mu_{i} は未知の独立な確率変数。
　　⇒現実的には、個別効果 \mu_{i} と説明変数 X_{i,t} が相関しないとは考えづらい。


ハウスマン検定

説明変数と誤差項の検定としてハウスマン検定がある。
帰無仮説が棄却されれば「固定効果モデル」、棄却できなければ「変量効果モデル」を採択する。問題点もあるようで今は使われない？

ハウスマン検定
　帰無仮説：固別効果と誤差項が独立である（相関しない）
　対立仮設：固別効果と誤差項が独立ではない（相関する）


プログラムコード

library(plm)
 
#-- 固定効果モデルの推定
fit1 <- plm(y ~ x1 + x2, data=ADS, model="within")
summary(fit1)
 
#-- 変量効果モデルの推定
fit2 <- plm(y ~ x1 + x2, data=ADS, model="random")
summary(fit2)
 
#--ハウスマン検定
phtest(fit1, fit2)

proc panel data=ADS ;
    id state date ;
    model y = x1 x2 s;
run;

参考

統計学（出版：東京図書）, 日本統計学会編
https://www.jil.go.jp/institute/zassi/backnumber/2015/04/pdf/006-009.pdf
SAS/ETS(R) 9.3 User's Guide

【SAS】Excelの文字切れ対策

SAS VBA

Excel文字切れ対策のめも　SASバージョン
VBAでの文字切れ対策記事(2020年11月)のあと直ぐ書きたかったが、3年以上寝かせていた。

基本的な考え方はVBAでの文字切れ対策記事を参照。
cochineal19.hatenablog.com

%let WB=ファイル名.xlsx;
%let WS=シート名;
%let FONTJ=ＭＳ ゴシック; /* 日本語フォント */
%let FONTE=Times New Roman; /* 英語フォント */
%let FONTSIZE=10; /* フォントサイズ */
%let STROW=;  /* 開始行 */
%let STCOL=;  /* 開始列 */
%let ENROW=;  /* 終了行 */
%let ENCOL=;  /* 終了列 */
   
*-- フォントサイズを１pt大きくして高さを自動調整する ;
filename exe dde "excel|system";
data _null_;
  file exe;
  put %unquote(%bquote('[workbook.activate("[&WB.]&WS.")]'));
  put %unquote(%bquote('[select("r&STROW.c&STCOL.:r&ENROW.c&ENCOL.")]'));
  put %unquote(%bquote('[font.properties("&FONTJ.",,%eval(&FONTSIZE. + 1))]'));
  put %unquote(%bquote('[font.properties("&FONTE.",,%eval(&FONTSIZE. + 1))]'));
  put %unquote(%bquote('[row.height(,"r&STROW.c&STCOL.:r&ENROW.c&ENCOL.",,3)]')); 
  put '[select("r1c1")]';
  put '[a1.r1c1(false)]';
run;
 
*-- get.cell(17)で行の高さを算出する ;
filename exe dde "excel|system";
data _null_;
  file exe;
  put '[workbook.next()]';
  put '[workbook.insert(3)]';
  put '[workbook.name(,"macro")]';
run;
 
filename xmacro1 dde "excel|macro!c1" notab;
data _null_;
  length CMD $2000.;
  file xmacro1 ;
  do i = 1 to &ENROW.;
    CMD = cats("=formula(get.cell(17,'&WS.", "'!rc), rc[1])");
    put CMD;
  end;
  put '=halt(true)';
  put '!dde_flush';
  file exe;
  put '[run("macro!r1c1:r1c1")]';
run;
 
*-- 行の高さを取り込む;
filename xmacro2 dde "excel|macro!r1c2:r%eval(&ENROW.+10)c2" notab;
data _HEIGHT1;
  infile xmacro2 ;
  length ROWC $20;
  ROW  = _n_;
  ROWC = cats("r",vvalue(ROW),"c1");
  input HEIGHT;
  HEIGHT2 = ceil(HEIGHT/15) * 15;
  if HEIGHT2 > 405 then HEIGHT2 = 405;

  if ROW >= &STROW.;
run;
 
*-- コード生成(DDE);
proc sort data=_HEIGHT1; by HEIGHT2 ROW; run;
data _HEIGHT1;
  set _HEIGHT1;
  by _HEIGHT2 ROW;
  if first.HEIGHT2 then ID=0; else ID+1;
  GRP = int(ID/10);
run;

proc transpose data=_HEIGHT1 out=_HEIGHT2;
  by HEIGHT2 GRP;
  var ROWC;
run;

data _HEIGHT3;
  set _HEIGHT2;
  length RNG1 DDE1 $20000;
  RNG1 = catx(",", of COL:);
  DDE1 = cats("put '[select(""", strip(RNG1), """)][row.height(", vvalue(HEIGHT2), ")]';");
run;

*-- DDE実行（高さ調整）;
proc sql noprint;
  select DDE1 into :DDE1 separated by ";" from _HEIGHT3;
quit;

filename exe dde "excel|system";
data _null_;
  file exe;
  put %unquote(%bquote('[workbook.activate("[&WB.]&WS.")]'));
  &DDE1.;
  put %unquote(%bquote('[select("r&STROW.c&STCOL.:r&ENROW.c&ENCOL.")]'));
  put %unquote(%bquote('[font.properties("&FONTJ.",,%eval(&FONTSIZE.))]'));
  put %unquote(%bquote('[font.properties("&FONTE.",,%eval(&FONTSIZE.))]'));
  put '[select("r1c1")]';

  %let ERR=ERR;
  put %unquote(%bquote('[&ERR.OR(false)]'));
  put %unquote(%bquote('[workbook.delete("macro")]'));
  put %unquote(%bquote('[&ERR.OR(true)]'));
run;

【統計】時系列分析（状態空間モデルをRで実行する、KFAS）

R、R Studio 統計

KFAS（カルマンフィルタ＆スムーザ）での状態空間モデルについてのメモ。その２

【目次】


以下で基本部分を取りまとめた、今回は実データで推定してみる。
【統計】時系列分析（状態空間モデル、KFAS） - こちにぃるの日記


使用するデータ

以前にARIMAモデルの記事で使用した横浜市の月別火災発生件数を題材とする。
【統計】時系列分析（ARIMA、SARIMA、SARIMAXをRで実行する） - こちにぃるの日記

ソース：
横浜市のオープンデータ
気象庁ホームページ

2000年1月～2019年12月の20年間240カ月分を用いる。
ARIMAモデルの記事と同様に、17年分（2000～2016年）を学習期間、３年（2017～2019年）をテスト期間する（ホールドアウト法）。
また、火災と関連がありそうなものとして、横浜市の月別平均気温を回帰成分（外生変数）として用いる。

#-- import
FIRE <- read_csv("./yokohama_kasai.csv")
TEMP <- read.csv("./yokohama_kion.csv")

#-- データ加工（時系列、回帰成分）
colnames(FIRE) <- c("YEARJ","YEAR","MONTH","FIRE")
ADTS <- data.frame(FIRE, TEMP=TEMP$平均気温)
ADTS$TESTFL <- ifelse(ADTS$YEAR >= 2017, 1, 0)  #-- 2017年以降をテストデータ
ADTS$CNT    <- ifelse(ADTS$TESTFL == 0, ADTS$FIRE, NA)  #-- テストデータはNA


状態空間モデルの実装

今回は、ローカル線形トレンドモデルを使って以下の３モデルを比較してみる。
・ローカル線形トレンドモデル
・ローカル線形トレンドモデル＋季節成分
・ローカル線形トレンドモデル＋季節成分＋回帰成分（平均気温）

#-- ローカル線形トレンド
mod1 <- SSModel(ADTS$CNT ~ SSMtrend(2, Q=list(NA,NA)), H=NA)
fit1 <- fitSSM(mod1, inits=numeric(3), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean","disturbance"))

#-- ローカル線形トレンド + 季節成分
mod2 <- SSModel(ADTS$CNT ~ SSMtrend(2, Q=list(NA,NA)) + SSMseasonal(12, Q=NA), H=NA)
fit2 <- fitSSM(mod2, inits=numeric(4), method="BFGS")
kfs2 <- KFS(fit2$model, smoothing=c("state","mean","disturbance"))

#-- ローカル線形トレンド + 季節成分 + 回帰成分
mod3 <- SSModel(ADTS$CNT ~ SSMtrend(2, Q=list(NA,NA)) + SSMseasonal(12, Q=NA) + SSMregression(~ ADTS$TEMP, Q=NA), H=NA)
fit3 <- fitSSM(mod3, inits=numeric(5), method="BFGS")
kfs3 <- KFS(fit3$model, smoothing=c("state","mean","disturbance"))


実行結果は以下の通り（学習データ：2000～2016年、予測：2017～2019年）。



モデル評価

残差分析は今回割愛する。本来は残差を見て外生変数の取捨選択や経時的な構造変化の有無などを検討する。
例えば、大きな災害や技術革新によって構造変化があれば、変化時点＝１、それ以外＝０としたダミー変数を設けるなどしてモデルの改善を目指す。

以下はテストデータに対するAICとRMSE、MAE、MAPEによる評価を記載。
長期予測の評価においては、AICは不向きなようである（AICは１期先）。
RMSE、MAE、MAPEなどの評価指標については、データの特性に応じて用いる指標を決める。
例えばスパイクの多いデータではRMSEは不向きかもしれない（スパイクに過剰に反応してしまうかもしれない）。0件が発生するデータではMAPEは不向きかもしれない（分母ゼロで計算できない）。

> #-- AIC
> logLik.1 <- kfs1$logLik - sum(kfs1$Finf>0) * log(2*pi)/2
> logLik.2 <- kfs2$logLik - sum(kfs2$Finf>0) * log(2*pi)/2
> logLik.3 <- kfs3$logLik - sum(kfs3$Finf>0) * log(2*pi)/2
> AIC1 <- -2*logLik.1 + 2*(3 + 1)           #-- 攪乱項(NA数)+トレンド成分
> AIC2 <- -2*logLik.2 + 2*(4 + 1 + 11)      #-- 攪乱項(NA数)+トレンド成分+季節成分
> AIC3 <- -2*logLik.3 + 2*(5 + 1 + 11 + 1)  #-- 攪乱項(NA数)+トレンド成分+季節成分+回帰成分 ←回帰成分の数え方はこれで良い？
> print(paste0(round(c(AIC1,AIC2,AIC3),4)))
[1] "1787.9992" "1640.0171" "1643.6485"

> #-- RMSE, MAE, MAPE
> ADTS2 <- ADTS %>% mutate(KFAS1=c(kfs1$muhat), KFAS2=c(kfs2$muhat), KFAS3=c(kfs3$muhat))
> ADTS3 <- ADTS2[ADTS2$TESTFL==1,] #-- テストデータのみにする
> 
> fnc1 <- function(KFAS){
+   RMSE <- sqrt(1/length(ADTS3$FIRE) * sum(KFAS - ADTS3$FIRE)^2)
+   MAE  <- 1/length(ADTS3$FIRE) * sum(abs(KFAS - ADTS3$FIRE))
+   MAPE <- 1/length(ADTS3$FIRE) * sum(abs(KFAS - ADTS3$FIRE) / ADTS3$FIRE) * 100
+   print(paste0(round(c(RMSE,MAE,MAPE),4)))
+ }
> fnc1(ADTS3$KFAS1)
[1] "133.343" "23.4793" "43.7605"
> fnc1(ADTS3$KFAS2)
[1] "2.1789"  "10.7435" "16.9743"
> fnc1(ADTS3$KFAS3)
[1] "0.0949"  "10.7321" "16.8962"


参考

https://cran.r-project.org/web/packages/KFAS/KFAS.pdf
カルマンフィルタ Rを使った時系列予測と状態空間モデル（出版：共立出版）, 野村俊一
https://www.actuaries.jp/lib/y_ronbun/pdf/2017.pdf
https://www.jstage.jst.go.jp/article/seitai/66/2/66_361/_pdf

【統計】時系列分析（状態空間モデル、KFAS）

統計 R、R Studio

KFAS（カルマンフィルタ＆スムーザ）での状態空間モデルについてのメモ。
まとめておきたかったが、ずっと手を付けられず。。
分かっていない部分も多いのであくまで個人備忘録として。
個人的に野村先生の書籍が大変参考になった。

【目次】


ARIMA関連は以下の記事
【統計】時系列分析（AR、MA、ARMA、ARIMA、SARIMA、ARIMAXモデル） - こちにぃるの日記
【統計】時系列分析（ARIMA、SARIMA、SARIMAXをRで実行する） - こちにぃるの日記
【統計】時系列分析（ARIMA、SARIMA、SARIMAXをSASで実行する） - こちにぃるの日記


状態空間モデル（State Space Model）

時系列データを解析することができる統計手法。
時点間の相関をデザインすることができ、将来予測にも用いることができる。

状態空間モデルでは、実際に得られる観測値と本来の状態を、それぞれ『観測モデル（観測方程式）』と『状態モデル（状態方程式）』に分けてモデリングする。
以下概念図の通り、本来の状態から何かしらのノイズが加わって観測値が得られる。




なお、推定目的によって、以下の用語の使い分けがあるとのこと。
・フィルタ（filtering）：現在の状態を推定する。
・平滑化（smoothing）：過去の状態を推定する。欠損値も推定。
・予測（forecasting）：将来の状態を推定する。


KFASでの実装

KFASでの実装は、Modeling ⇒ Fitting ⇒ Filtering&Smoothing の順で行い、それぞれ SSModel，fitSSM，KFS関数を使用する。
以下はローカルレベルモデルの実装例。SSMtrend関数でトレンド成分を設定する。

library(KFAS)
mod1 <- SSModel(ADS$Y ~ SSMtrend(1, Q=NA), H=NA))
fit1 <- fitSSM(mod1, inits=numeric(2), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean"))


季節成分や回帰成分（外生変数）を含める場合は、SSMseasonal関数と SSMregression関数で設定する。
また、AR成分もSSMarima関数で設定できるようだが、SSMarimaは使用したことはない。

library(KFAS)
mod1 <- SSModel(ADS$Y ~ SSMtrend(1, Q=NA)
                          + SSMseasonal(12, Q=NA)
                          + SSMregression(~ ADS$VAR1 + ADS$VAR2, Q=NA)
                          ,H=NA)
fit1 <- fitSSM(mod1, inits=numeric(4), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean"))


KFASの特徴は、ポアソン分布や二項分布などの正規分布（ガウシアン）以外も使用できることにある。
以下はポアソン分布を指定した例。シミュレーションが入るのでPCが弱々だと厳しい・・・

mod1p <- SSModel(ADTS$CNT ~ SSMtrend(1, Q=NA)
                 + SSMseasonal(12, Q=NA)
                 + SSMregression(~ ADTS$TEMP, Q=NA)
                 ,distribution="poisson"
                 ,u=1)
fit1p <- fitSSM(mod1p, inits=numeric(3), method="BFGS", nsim=1000)
kfs1p <- KFS(fit1p$model, nsim=1000)


状態空間モデルで扱われるトレンド成分モデル

■ ローカルレベルモデル

上記概念図で示した期待値がランダムウォークすると仮定したモデル。
「１時点前の期待値＋ホワイトノイズ」。
\quad y _{t}= \mu _{t} + \varepsilon _{t} \ , \ \varepsilon_{t}\sim 𝑁(0, H)
\quad \mu _{t}=\mu _{t-1}+\eta _{t-1} \ , \ \eta _{t-1}\sim 𝑁(0, 𝑄) 

mod1 <- SSModel(ADS$Y ~ SSMtrend(1, Q=NA), H=NA)
fit1 <- fitSSM(mod1, inits=numeric(2), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean"))



■ 平滑化トレンドモデル

傾きがランダムウォークすると仮定したモデル。
「１時点前と２時点前の階差系列＋ホワイトノイズ」。
\quad y _{t}= \mu _{t} + \varepsilon _{t} \ , \ \varepsilon_{t}\sim 𝑁(0, H)
\quad \mu _{t}=\mu _{t-1} + (\Delta \mu _{t-1} + \eta _{t-1}) \ , \ \eta _{t-1}\sim 𝑁(0, 𝑄)
\quad \Delta \mu _{t-1} = \mu _{t-1} - \mu _{t-2} 

mod1 <- SSModel(ADS$Y ~ SSMtrend(2, Q=list(0,NA)), H=NA)
fit1 <- fitSSM(mod1, inits=numeric(2), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean"))



■ ローカル線形トレンドモデル

期待値と傾き両方がランダムウォークすると仮定したモデル。
\quad y _{t}= \mu _{t} + \varepsilon _{t} \ , \ \varepsilon_{t}\sim 𝑁(0, H)
\quad \mu _{t}=(\mu _{t-1} + \eta _{t-1,1}) + (\Delta \mu _{t-1} + \eta _{t-1,2}) \ , \ \eta _{t-1,1}\sim 𝑁(0, 𝑄1) \ , \ \eta _{t-1,2}\sim 𝑁(0, 𝑄2)
\quad \Delta \mu _{t-1} = \mu _{t-1} - \mu _{t-2} 

mod1 <- SSModel(ADS$Y ~ SSMtrend(2, Q=list(NA,NA)), H=NA)
fit1 <- fitSSM(mod1, inits=numeric(3), method="BFGS")
kfs1 <- KFS(fit1$model, smoothing=c("state","mean"))


参考

https://cran.r-project.org/web/packages/KFAS/KFAS.pdf
カルマンフィルタ Rを使った時系列予測と状態空間モデル（出版：共立出版）, 野村俊一
https://www.actuaries.jp/lib/y_ronbun/pdf/2017.pdf
https://www.jstage.jst.go.jp/article/seitai/66/2/66_361/_pdf

【R言語＆Batch】Rスクリプトをバッチ実行する

R、R Studio コマンドプロンプト

これまでJscriptを使ってRを実行する方法は以下記事で触れていたが、バッチ(.bat)での実行方法の記事は無かったので作成。
cochineal19.hatenablog.com
cochineal19.hatenablog.com

Rスクリプトをこのバッチファイルにドラッグ＆ドロップするとRが実行され、ログファイルが出力される仕様。
下記コード内のRscript_exeは各個人のexe保存先を設定。
R側の受け取り方法は上記その２を参照。

@echo off

rem:===================================================================
rem: Setting
rem:===================================================================

set Rscript_exe="C:\Program Files\R\R-4.2.3\bin\Rscript.exe"

set RscriptPath="%~dp1%~nx1"
set RscriptPath=%RscriptPath:\=/%

set RscriptFolder="%~dp1"
set RscriptFolder=%RscriptFolder:\=/%

set RscriptName="%~nx1"

set yymmdd=%date:/=%
set yymmdd=%yymmdd: =0%

set LogName="%~dp1__%~n1_%yymmdd%.log"
echo %LogName%


rem:===================================================================
rem: Exec
rem:===================================================================

echo 以下がドラッグ＆ドロップされました。
echo ・フォルダ名：%RscriptFolder:"=%
echo ・ファイル名：%RscriptName:"=%
echo ・拡張子名　：%~x1
echo.

set TRUE_FALSE=FALSE
if "%~x1"==".R" set TRUE_FALSE=TRUE
if "%~x1"==".r" set TRUE_FALSE=TRUE

if %TRUE_FALSE%==FALSE (
  echo Rスクリプトをドラッグ＆ドロップしてください。
  goto no1
)


set /P selected="実行しますか？ 実行する場合、1を入力してください。"
if /i %selected%==1 (
  goto yes1
)


:no1
echo 処理を中止しました。
goto exitlabel


:yes1
echo 実行します。


echo ------------------------------------------------------------- > %LogName%
echo 　対象ファイル：%RscriptPath% >> %LogName%
echo 　実行日時　　：%RscriptPath% >> %LogName%
echo ------------------------------------------------------------- >> %LogName%

%Rscript_exe% --slave --vanilla %RscriptPath% %RscriptFolder% "__error.txt" >> %LogName% 2>&1

echo -- >> %LogName%

echo ------------------------------------------------------------- >> %LogName%
echo 　End of File >> %LogName%
echo ------------------------------------------------------------- >> %LogName%


echo 終了しました。

:exitlabel
pause

rem:===================================================================
rem: End of File
rem:===================================================================

【統計】時系列分析（ARIMA、SARIMA、SARIMAXをSASで実行する）

SAS 統計

ARIMAプロシジャのメモ。
基本はRを使っているので初歩。

【目次】


以前の記事
【統計】時系列分析（AR、MA、ARMA、ARIMA、SARIMA、ARIMAXモデル） - こちにぃるの日記
【統計】時系列分析（ARIMA、SARIMA、SARIMAXをRで実行する） - こちにぃるの日記


Auto ARIMA

identify ステートメントでAR, MA過程の複数パターンを評価できる。
var=target で目的変数を指定。カッコ内に階差を指定。
nlag で季節サイクルを指定。
p=(0:10) q=(0:10) でAR過程を０～10の各組合せで推定。
stationarity=(adf) で定常性を検定。adf で拡張ディッキー–フラー検定。
crosscorr=(variable_name) で外生変数（回帰成分）との相関を評価。
esacf minic scan あたりで各組合せを評価。

proc arima data=adts ;
    identify var=target(1,12) nlag=12
             p=(0:10) q=(0:10)
             stationarity=(adf)
             crosscorr=(variable_name) outcov=_outcov
             esacf minic scan ;
run;
quit;

ARIMA(p,0,0)：AR(p)モデル

estimate ステートメントでAR過程とMA過程を指定。

proc arima data=adts ;
    identify var=target;
    estimate p=2; 
run;
quit;

ARIMA(0,0,q)：MA(q)モデル

proc arima data=adts ;
    identify var=target;
    estimate q=2; 
run;
quit;

ARIMA(p,0,q)：ARMA(p,q)モデル

proc arima data=adts ;
    identify var=target;
    estimate p=2 q=2; 
run;
quit;

ARIMA(p,d,q)：ARIMA(p,d,q)モデル

proc arima data=adts ;
    identify var=target(1);
    estimate p=2 q=2; 
run;
quit;

ARIMAXモデル

identify ステートメントの crosscorr と estimate ステートメントの input で外生変数を指定。

proc arima data=adts ;
    identify var=target(1) crosscorr=(variable_name);
    estimate p=2 q=2 input=(variable_name);
run;
quit;

SARIMAモデル

nlag でサイクルを指定。

proc arima data=adts ;
    identify var=target(1,12) nlag=12;
    estimate p=(1 2)(1 2) q=(1 2)(1 2); 
run;
quit;

SARIMAXモデル

proc arima data=adts ;
    identify var=target(1,12) nlag=12 crosscorr=(variable_name);
    estimate p=(1 2)(1 2) q=(1 2)(1 2) input=(variable_name);
run;
quit;

将来予測

forecast ステートメントで将来予測。
lead で何期まで予測するか指定。
interval でサイクルの単位。
id で inputデータの必要な変数を outputデータに保持させる。

proc arima data=adts ;
    identify var=target(1,12) nlag=12 crosscorr=(variable_name);
    estimate p=(1 2)(1 2) q=(1 2)(1 2) input=(variable_name);
    forecast lead=12 interval=month id=id_variable_name out=_pred ;
run;
quit;
  
proc sgplot data=_pred;
   band Upper=u95 Lower=l95 x=id_variable_name / LegendLabel="95%CI";
   scatter x=id_variable_name y=target;
   series  x=id_variable_name y=forecast ;
run;


参考

https://support.sas.com/documentation/onlinedoc/ets/132/arima.pdf

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。