TSA007 - waseda2bun_r @ ウィキ

#実質・原系列 x <- read.csv("gdp_ex.csv",skip=3) gdp <- x[,2] ex <- x[,3] GDP <- ts(log(gdp),start=c(1980,1),frequency=4)#log EX <- ts(log(ex),start=c(1980,1),frequency=4)#log arima_g1 <- arima(EX,order=c(2,1,2),transform.pars =FALSE);arima_g1 arima_g2 <- arima(EX,order=c(2,1,2));arima_g2 #arima_g2はARの係数が定常性を満たす範囲内で尤度を最大化する arima_ex1 <- arima(EX,order=c(2,1,2),transform.pars =FALSE);arima_ex1 arima_ex2 <- arima(EX,order=c(2,1,2));arima_ex2 #arima_ex2はARの係数が定常性を満たす範囲内で尤度を最大化する

 ・結果

Call: arima(x = GDP, order = c(2, 1, 2), transform.pars = FALSE) Coefficients:          ar1      ar2      ma1     ma2       0.0055  -0.9344  -0.5808  0.8827 s.e.  0.0355   0.0334   0.0628  0.0392 sigma^2 estimated as 0.001039:  log likelihood = 241.34,  aic = -472.68       Call: arima(x = GDP, order = c(2, 1, 2)) Coefficients:          ar1      ar2      ma1     ma2       0.0055  -0.9345  -0.5808  0.8827 s.e.  0.0355   0.0335   0.0628  0.0392 sigma^2 estimated as 0.001039:  log likelihood = 241.34,  aic = -472.68

警告メッセージ： In log(s2) :  計算結果が NaN になりました Call: arima(x = EX, order = c(2, 1, 2), transform.pars = FALSE) Coefficients:          ar1      ar2      ma1     ma2       0.0602  -0.9811  -0.1803  1.0002 s.e.  0.0225   0.0159   0.0221  0.0578 sigma^2 estimated as 0.002589:  log likelihood = 185.48,  aic = -360.95     Call: arima(x = EX, order = c(2, 1, 2)) Coefficients:          ar1      ar2     ma1     ma2       0.0118  -0.9996  0.0062  0.9993 s.e.  0.0085   0.0012  0.0222  0.0517 sigma^2 estimated as 0.002479:  log likelihood = 188.2,  aic = -366.4

 ・季節調整をする

<- read.csv("gdp_ex.csv",skip=3) gdp <- x[,2] ex <- x[,3] GDP <- ts(log(gdp),start=c(1980,1),frequency=4)#log EX <- ts(log(ex),start=c(1980,1),frequency=4)#log arima_g1 <- arima(GDP,order=c(2,1,2),transform.pars =FALSE);arima_g1 arima_g2 <- arima(GDP,order=c(2,1,2));arima_g2 #arima_g2はARの係数が定常性を満たす範囲内で尤度を最大化する arima_ex1 <- arima(EX,order=c(2,1,2),transform.pars =FALSE);arima_ex1 arima_ex2 <- arima(EX,order=c(2,1,2));arima_ex2 #arima_ex2はARの係数が定常性を満たす範囲内で尤度を最大化する sa <- list(order=c(1,1,1),period=4) sarima_g1 <- arima(GDP,order=c(2,1,2),seasonal=sa,transform.pars =FALSE);arima_g1 sarima_g2 <- arima(GDP,order=c(2,1,2),seasonal=sa);arima_g2 sarima_e1 <- arima(EX,order=c(2,1,2),seasonal=sa,transform.pars =FALSE) sarima_e2 <- arima(EX,order=c(2,1,2),seasonal=sa) sarima_e1 sarima_e2 par(mfrow=c(2,1)) ahat <- sarima_e1$resid EX_hat <- EX-ahat plot(EX,type="l",main="log EXO") lines(EX_hat,lty=2,col=2) plot(ahat,type="l",col=4,main="Resid") abline(h=mean(ahat)) lines(EX_hat,lty=2,col=2)

・SARIMAモデルの検定

b1 <- arima_sa_e1$coef V1 <- arima_sa_e1$var.coef t1 <- numeric(4) for (j in 1:4){ t1[j] <- b[j]/sqrt(V[j,j])} names(t) <- c("t_ar1","t_ar2","t_ma2","t_ma2");t test1 <- ((t1<0)&(pnorm(t1)<0.05))|((t1>0)&(pnorm(t1)>0.95));test1 b2 <- arima_sa_g1$coef V2 <- arima_sa_g1$var.coef t2 <- numeric(4) for (j in 1:4){ t2[j] <- b[j]/sqrt(V[j,j])} names(t2) <- c("t_ar1","t_ar2","t_ma2","t_ma2");t test2 <- ((t2<0)&(pnorm(t2)<0.05))|((t2>0)&(pnorm(t2)>0.95));test2 st.test1 <- abs(polyroot (c(-b1[2],-b1[1],1)));st.test1 st.test2_1 <- (st.test1<1);st.test2_1 st.test2 <- abs(polyroot (c(-V1[2],-V1[1],1)));st.test2 st.test2_2 <- (st.test2<1);st.test2_2

・残差の検定

ahat <- sarima_e1$resid Box.test(ahat,lag=1,type="Ljung-Box") #h0=残差系列に自己相関がない。 #0.05以上→5％で棄却されない＝ホワイトノイズとみなせる ahat1 <- sarima_e1$resid Box.test(ahat1,lag=1,type="Ljung-Box") #h0=残差系列に自己相関がない ahat2 <- sarima_g2$resid Box.test(ahat2,lag=1,type="Ljung-Box") #h0=残差系列に自己相関がない。 #0.05以上→5％で棄却されない＝ホワイトノイズとみなせる ahat <- sarima_e1$resid Box.test(ahat,lag=1,type="Ljung-Box") #h0=残差系列に自己相関がない

・検定

pval <- numeric(20) for (j in 1:20){ kekka <- Box.test(ahat1,lag=j,type="Ljung-Box") pval[j] <- kekka$p.value } ttl <- "Ljung-Box test p value"; y1 <- c(0,1) plot(pval,xlab="s",ylim=c(0.05,1),main=ttl) abline(h=0.05,col=2) #保存されない tsdiag(sarima_g2);tsdiag(sarima_g2,lag=20) tsdiag(sarima_e2);tsdiag(sarima_e2,lag=20)