VAR analysis 2017Q3

kikan <- c("1992::2017-09-30")
getSymbols("RSXFS",src="FRED")
getSymbols("RRSFS",src="FRED")
getSymbols("PERMIT",src="FRED")
getSymbols("HOUST",src="FRED")
beep(2)
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
getSymbols("GDPC1",src="FRED")
getSymbols("QUSR628BIS",src="FRED");last(QUSR628BIS,n=3)
PROP <- append(QUSR628BIS,as.xts(c(125.84,127.84),c(as.Date("2017-04-01"),as.Date("2017-07-01"))))


v_GPC_q_1992_2017 <- c()
v_GPC_q_1992_2017 <- merge(GDPC1[kikan],to.quarterly(UNDCONTSA[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(PAYEMS[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(RRSFS[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(RSXFS[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(PERMIT[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(HOUST[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(SP5[kikan])[,4])
# v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(QUSR628BIS[kikan])[,4])
v_GPC_q_1992_2017 <- merge(v_GPC_q_1992_2017,to.quarterly(PROP[kikan])[,4])


names(v_GPC_q_1992_2017)[2] <- "UNDCONTSA"
names(v_GPC_q_1992_2017)[3] <- "PAYEMS"
names(v_GPC_q_1992_2017)[4] <- "RRSFS"
names(v_GPC_q_1992_2017)[5] <- "RXSFS"
names(v_GPC_q_1992_2017)[6] <- "PERMIT"
names(v_GPC_q_1992_2017)[7] <- "HOUST"
names(v_GPC_q_1992_2017)[8] <- "SP5"
names(v_GPC_q_1992_2017)[9] <- "PROP"

VARselect(v_GPC_q_1992_2017)
l <- 16
m_aic <- VARselect(v_GPC_q_1992_2017)$selection[1]
predict(VAR(v_GPC_q_1992_2017,la=m_aic),n.ahead = 16)$fcst
plot(seq(as.Date(strsplit(kikan,"::")[[1]][2]),as.Date("2100-10-01"),by="quarters")[2:(l+1)],predict(VAR(v_GPC_q_1992_2017,la=m_aic),n.ahead = l)$fcst$PAYEMS[,1])

do forecast by ARIMA

The length of periods are given as the parameter "h".


library(forecast)
am <- auto.arima( PA[k1992], ic="aic", trace=T, stepwise=F, approximation=F, start.p=0, start.q=0, start.P=0, start.Q=0)
forecast(am,level = c(50,95), h = 12)

VAR - 1992::2016-12-31

prepare data for extended period.

kikan <- c("1992::2016-12-31")
v_GPC_q_1992_2016 <- merge(GDPC96[kikan],to.quarterly(UNDCONTSA[kikan])[,4])
v_GPC_q_1992_2016 <- merge(v_GPC_q_1992_2016,to.quarterly(PAYEMS[kikan])[,4])
names(v_GPC_q_1992_2016)[2] <- "UNDCONTSA"
names(v_GPC_q_1992_2016)[3] <- "PAYEMS"

> VARselect(v_GPC_q_1992_2016)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
    10      2      2     10 

$criteria
                  1            2            3            4            5            6            7            8            9           10
AIC(n) 2.785261e+01 2.545883e+01 2.543818e+01 2.542177e+01 2.544280e+01 2.546334e+01 2.552722e+01 2.558417e+01 2.532732e+01 2.523473e+01
HQ(n)  2.798702e+01 2.569405e+01 2.577421e+01 2.585860e+01 2.598043e+01 2.610179e+01 2.626647e+01 2.642423e+01 2.626818e+01 2.627640e+01
SC(n)  2.818592e+01 2.604212e+01 2.627145e+01 2.650502e+01 2.677603e+01 2.704655e+01 2.736041e+01 2.766735e+01 2.766047e+01 2.781786e+01

FPE(n) 1.248282e+12 1.140354e+11 1.119086e+11 1.104548e+11 1.134074e+11 1.166774e+11 1.257451e+11 1.350810e+11 1.064882e+11 9.943996e+10

This will give the predict from now to 16 quarter latter.

predict(VAR(v_GPC_q_1992_2016,la=10),n.ahead = 16)

New GDP VAR model based on quarterly data.

New S&P 500 forecast is based upon quarterly data. Therefore new quarterly data based VAR model for GDP is required.

v_GPC_q_1992_2016 <- merge(GDPC96["1992::2016-04-01"],to.quarterly(UNDCONTSA["1992::2016-04-01"])[,4])
v_GPC_q_1992_2016 <- merge(v_GPC_q_1992_2016,to.quarterly(PAYEMS["1992::2016-04-01"])[,4])
names(v_GPC_q_1992_2016)[2] <- "UNDCONTSA"
names(v_GPC_q_1992_2016)[3] <- "PAYEMS"
VARselect(v_GPC_q_1992_2016)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
    10      2      2     10 

$criteria
                  1            2            3            4            5            6            7            8            9
AIC(n) 2.789256e+01 2.551374e+01 2.550815e+01 2.549235e+01 2.552063e+01 2.551927e+01 2.558270e+01 2.561892e+01 2.537707e+01
HQ(n)  2.802866e+01 2.575191e+01 2.584840e+01 2.593467e+01 2.606502e+01 2.616573e+01 2.633124e+01 2.646954e+01 2.632976e+01
SC(n)  2.823038e+01 2.610492e+01 2.635270e+01 2.659026e+01 2.687190e+01 2.712390e+01 2.744070e+01 2.773029e+01 2.774180e+01
FPE(n) 1.299168e+12 1.204798e+11 1.200424e+11 1.185840e+11 1.226860e+11 1.235582e+11 1.332081e+11 1.403119e+11 1.124435e+11
                 10
AIC(n) 2.525669e+01
HQ(n)  2.631145e+01
SC(n)  2.787478e+01
FPE(n) 1.023095e+11

Then "lag.max" parameter for the new model is 10.

predict(VAR(v_GPC_q_1992_2016,lag.max=10))
my_sp5(predict(VAR(v_GPC_q_1992_2016,lag.max=10))$fcst$GDPC96[,1],predict(VAR(v_GPC_q_1992_2016,lag.max=10))$fcst$PAYEMS[,1],predict(VAR(v_GPC_q_1992_2016,lag.max=10))$fcst$UNDCONTSA[,1])

private function "my_sp5" will be explained in the  future post .

GDP forecast by VAR

Thus we should be able to see updated GDP forecast. Most of the figures are almost not changed at all.

> predict(VAR(v_XCPG_1992_2016),lag=6)$fcst$GDP
          fcst    lower    upper        CI
 [1,] 16633.79 16582.69 16684.89  51.09898
 [2,] 16653.44 16581.13 16725.76  72.31579
 [3,] 16672.98 16584.34 16761.61  88.63489
 [4,] 16692.39 16589.96 16794.82 102.42740
 [5,] 16711.69 16597.08 16826.30 114.61098
 [6,] 16730.88 16605.22 16856.53 125.65568
 [7,] 16749.97 16614.13 16885.81 135.84034
 [8,] 16768.95 16623.61 16914.30 145.34629
 [9,] 16787.85 16633.55 16942.15 154.29927
[10,] 16806.65 16643.86 16969.44 162.79063

Updated dataframe for VAR analysis

use GDP,RSXFS, UNDCONTSA and PAYEMS.

v_XCPG_1992_2016 <- merge(RSXFS["1992::2016-04-01"],UNDCONTSA["1992::2016-04-01"])
v_XCPG_1992_2016 <- merge(v_XCPG_1992_2016,PAYEMS["1992::2016-04-01"])
v_XCPG_1992_2016 <- merge(v_XCPG_1992_2016,GDP_XTS["1992::2016-04-01"])

change column name.

> tail(v_XCPG_1992_2016)
            RSXFS UNDCONTSA PAYEMS GDP_XTS..1992..2016.04.01..
2015-11-01 395480       964 142875                    16485.24
2015-12-01 396723       976 143146                    16499.92
2016-01-01 394878       976 143314                    16514.59
2016-02-01 395126       987 143547                    16547.33
2016-03-01 393775       991 143733                    16580.67
2016-04-01 399134      1001 143877                    16614.00
> names(v_XCPG_1992_2016)[4]
[1] "GDP_XTS..1992..2016.04.01.."
> names(v_XCPG_1992_2016)[4] <- "GDP"

plot past 60 months and next 10 months.

>plot(c(as.numeric(last(GDP_XTS,n=60)),predict(VAR(v_XCPG_1992_2016),lag=6)$fcst$GDP[,1]),type="l")

re-re-forecast GDP

Based upon new update GDPnow forecast at May 28, do re-re-forecast GDP for next 10 months. 2nd qtr. New GDP annual growth rate is 2.8%.

> predict(VAR(v_XCPG_1992_2016,lag.max = 6))$fcst$GDP
          fcst    lower    upper        CI
 [1,] 16650.59 16620.69 16680.49  29.90123
 [2,] 16678.36 16613.26 16743.45  65.09441
 [3,] 16701.88 16596.66 16807.10 105.22027
 [4,] 16726.01 16586.58 16865.44 139.42727
 [5,] 16746.25 16575.42 16917.08 170.83343
 [6,] 16764.87 16565.62 16964.12 199.25079
 [7,] 16782.19 16554.59 17009.78 227.59567
 [8,] 16799.61 16544.35 17054.86 255.25926
 [9,] 16816.25 16533.71 17098.78 282.53532
[10,] 16832.42 16523.43 17141.40 308.98703

Replace a specific column in the data structure.

Prepare updated data in "w". v_XCPG_1992_2016 only has data from 1992. Use "w["1992::"]"  to limit data just from 1992. Data were updated according to GDPnow 2016/5/28.


> tail(v_XCPG_1992_2016$GDP)
2015-11-01 16477.960
2015-12-01 16485.354
2016-01-01 16505.700
2016-02-01 16543.000
2016-03-01 16576.000
2016-04-01 16607.000
> tail(w)
2015-11-01 16482.07
2015-12-01 16493.57
2016-01-01 16505.07
2016-02-01 16540.61
2016-03-01 16580.61
2016-04-01 16620.61
> v_XCPG_1992_2016$GDP <- w["1992::"]
> tail(v_XCPG_1992_2016$GDP)
                GDP
2015-11-01 16482.07
2015-12-01 16493.57
2016-01-01 16505.07
2016-02-01 16540.61
2016-03-01 16580.61
2016-04-01 16620.61

S&P500 forecast next 10 months

With newly updated GDP data and model at  May 25th , S&P500 forecast next 10 months are as below.

Please note that the current price is already more than 100 pts below theoretical value. 2323.292 is expected  at the end of this year.


> 0.603*(predict(VAR(v_XCPG_1992_2016,lag=6))$fcst$GDP[,1])-7791
 [1] 2238.212 2254.216 2268.441 2282.354 2293.338 2303.527 2313.299 2323.292 2332.697 2341.712

GDP forecast by newly updated data.

2016 1st Quarter GDP was revised upward. The growth rate was updated from 0.5% to 0.8%.
The new GDP forecast are as below. The previous one is available at May 23rd entry. #10 months is up 6.1 trillion #8 is also up 6.7. It will also give S&P 500 from another 30 to 40 pts at the end of this year.



> predict(VAR(v_XCPG_1992_2016,lag=6))$fcst$GDP
          fcst    lower    upper        CI
 [1,] 16632.19 16602.29 16662.10  29.90217
 [2,] 16658.73 16593.67 16723.80  65.06516
 [3,] 16682.32 16577.21 16787.44 105.11452
 [4,] 16705.40 16566.14 16844.65 139.25646
 [5,] 16723.61 16552.99 16894.24 170.62578
 [6,] 16740.51 16541.47 16939.54 199.03576
 [7,] 16756.71 16529.33 16984.10 227.38612
 [8,] 16773.29 16518.24 17028.34 255.05082
 [9,] 16788.88 16506.56 17071.21 282.32406
[10,] 16803.83 16495.07 17112.60 308.76632

predict by VAR analysis

In some cases, as VAR predict function's output  all implied data on the console, it may take extra time to find the data of focus. Specifying attribute as below suppresses unnecessary outputs.

>predict(VAR(v_XCPG_1992_2016),lag=6)$fcst$<tag of focus data>

Below will help to understand output's data structure.

>predict(VAR(v_XCPG_1992_2016),lag=6)[1]
>predict(VAR(v_XCPG_1992_2016),lag=6)[2]
>predict(VAR(v_XCPG_1992_2016),lag=6)[3]

Please see the sample outputs.

> predict(VAR(v_XCPG_1992_2016),lag=6)$fcst$GDP
          fcst         lower         upper        CI
 [1,] 16577.50 16526.29 16628.71  51.21148
 [2,] 16596.18 16523.71 16668.65  72.46879
 [3,] 16614.76 16525.95 16703.58  88.81411
 [4,] 16633.27 16530.64 16735.89 102.62429
 [5,] 16651.68 16536.86 16766.50 114.81930
 [6,] 16670.02 16544.15 16795.89 125.87053
 [7,] 16688.28 16552.22 16824.34 136.05765
 [8,] 16706.46 16560.90 16852.03 145.56265
 [9,] 16724.58 16570.07 16879.09 154.51178
[10,] 16742.63 16579.63 16905.63 162.99681
> 

how to do VAR analysis

まず、次元数を選択する。
> VARselect(v_RCPG_1992_2016)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     6      5      3      6 

$criteria
                  1            2            3            4            5            6            7
AIC(n) 3.559669e+01 3.416290e+01 3.382555e+01 3.367838e+01 3.358712e+01 3.355134e+01 3.356827e+01
HQ(n)  3.570054e+01 3.434984e+01 3.409558e+01 3.403150e+01 3.402332e+01 3.407063e+01 3.417064e+01
SC(n)  3.585565e+01 3.462902e+01 3.449884e+01 3.455884e+01 3.467474e+01 3.484613e+01 3.507023e+01
FPE(n) 2.880391e+15 6.867470e+14 4.901920e+14 4.232468e+14 3.865331e+14 3.732358e+14 3.800110e+14
                  8            9           10
AIC(n) 3.361131e+01 3.365361e+01 3.368294e+01
HQ(n)  3.429676e+01 3.442215e+01 3.453456e+01
SC(n)  3.532043e+01 3.556990e+01 3.580639e+01
FPE(n) 3.972810e+14 4.151925e+14 4.285075e+14

VAR分析を行う。VARselectで得られた最適次元数を指定する。ついでに予測を出力するためpredict関数を使用する。

> predict(VAR(v_RCPG_1992_2016,lag.max = 6))
$RRSFS
          fcst    lower    upper       CI
 [1,] 189476.9 186927.0 192026.9 2549.913
 [2,] 189763.4 186712.7 192814.2 3050.791
 [3,] 189747.8 186098.8 193396.9 3649.043
 [4,] 190150.1 185985.9 194314.2 4164.156

<出力結果中略>

$GDP
          fcst    lower    upper        CI
 [1,] 16592.21 16562.29 16622.13  29.91837
 [2,] 16616.42 16551.43 16681.41  64.99146
 [3,] 16642.72 16537.46 16747.98 105.25868
 [4,] 16667.69 16527.70 16807.69 139.99823
 [5,] 16695.66 16523.53 16867.80 172.13373
 [6,] 16722.72 16521.30 16924.13 201.41477
 [7,] 16749.10 16518.42 16979.77 230.67647
 [8,] 16774.03 16514.56 17033.50 259.47049
 [9,] 16798.65 16510.69 17086.60 287.95277
[10,] 16822.57 16507.10 17138.05 315.47539

> 

how to prepare structure for VAR analisys

Prepare data. each data should be XTS class.

> class(RRSFS)
[1] "xts" "zoo"

merge関数を使ってXTS classのデータを次々結合していく。

> v_RCPG_1992_2016 <- merge(RRSFS["1992::2016-03-01"],UNDCONTSA["1992::2016-03-01"])
> v_RCPG_1992_2016 <- merge(v_RCPG_1992_2016,PAYEMS["1992::2016-03-01"])
> v_RCPG_1992_2016 <- merge(v_RCPG_1992_2016,GDP_XTS["1992::2016-03-01"])
> colnames(v_RCPG_1992_2016)
[1] "RRSFS"                       "UNDCONTSA"                   "PAYEMS"                     
[4] "GDP_XTS..1992..2016.03.01.."
> tail(v_RCPG_1992_2016)
            RRSFS UNDCONTSA PAYEMS GDP_XTS..1992..2016.03.01..
2015-10-01 187883       943 142595                    16470.57
2015-11-01 188187       964 142875                    16477.96
2015-12-01 189231       976 143146                    16485.35
2016-01-01 188168       976 143314                    16492.75
2016-02-01 189012       987 143547                    16525.74
2016-03-01 188220       994 143755                    16558.72

四列目の列名が見にくいので変更する。

> names(v_RCPG_1992_2016)[4] <- c("GDP")
> head(v_RCPG_1992_2016)
            RRSFS UNDCONTSA PAYEMS      GDP
1992-01-01 118643       639 108377 9122.950
1992-02-01 118514       630 108314 9156.482
1992-03-01 117719       656 108368 9190.013
1992-04-01 118192       650 108527 9223.545
1992-05-01 118552       650 108653 9253.433

1992-06-01 118557       644 108721 9283.320