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
> 

update fundamental datas 統計 主要統計データ FRED YAHOO 指標

update data by getsymbols.

PAYEMS is All Employees: Total Nonfarm Payrolls
RRSFS is Real Retail and Food Services Sales
UNDCONTSA is New Privately-Owned Housing Units Under Construction: Total
PI is Personal Income
TWEXB is Trade Weighted U.S. Dollar Index: Broad

each operation will give data in XTS format.

getSymbols("GDPC96",src="FRED")
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED");
getSymbols("RRSFS",src="FRED")
getSymbols("PI",src="FRED")
getSymbols("TWEXB",src="FRED")

# PCE

getSymbols("PCEPI",src="FRED",auto.assign=TRUE)

# CPI

getSymbols("CPIAUCSL",src="FRED",auto.assign=TRUE)

# Core CPI

getSymbols("CPILFESL",src="FRED",auto.assign=TRUE)

# core PCE

getSymbols("PCEPILFE",src="FRED",auto.assign=TRUE)

how to convert from CSV to XTS

when time series data is given in CSV format as below.

~$ head SP5.csv 

Date,Open,High,Low,Close,Volume,Adj Close

2016-05-02,2067.169922,2084.870117,2025.910034,2052.320068,4068314600,2052.320068

2016-04-01,2056.620117,2111.050049,2033.800049,2065.300049,4087129000,2065.300049

2016-03-01,1937.089966,2072.209961,1937.089966,2059.73999,4379759000,2059.73999

the combination of read.csv, read.zoo and as.xts is able to transform and conver t data into xts objects.

> SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))

> tail(SP5)

              Open    High     Low   Close     Volume Adj.Close

2015-12-01 2082.93 2104.27 1993.26 2043.94 3922935900   2043.94

2016-01-04 2038.20 2038.20 1812.29 1940.24 5153017800   1940.24

2016-02-01 1936.94 1962.96 1810.10 1932.23 4881887000   1932.23

2016-03-01 1937.09 2072.21 1937.09 2059.74 4379759000   2059.74

2016-04-01 2056.62 2111.05 2033.80 2065.30 4087129000   2065.30

2016-05-02 2067.17 2084.87 2025.91 2052.32 4068314600   2052.32

> class(SP5)

[1] "xts" "zoo"

plot() function and frequency conversion function like to.monthly() are available as usual.

> plot(to.yearly(SP5))

access to attributes

xtsAttributes provides access to xts class specific attributes.

> xtsAttributes(PAYEMS)
$src
[1] "FRED"

$updated
[1] "2016-05-18 21:30:19 JST"

use attributes for base class.

> attributes(PAYEMS)
$dim
[1] 928   1

$index
  [1] -978307200 -975628800 -973209600 -970531200 -967939200 -965260800 -962668800 -959990400 -957312000
 [10] -954720000 -952041600 -949449600 -946771200 -944092800 -941587200 -938908800 -936316800 -933638400
 [19] -931046400 -928368000 -925689600 -923097600 -920419200 -917827200 -915148800 -912470400 -910051200

<skip>

[910] 1412121600 1414800000 1417392000 1420070400 1422748800 1425168000 1427846400 1430438400 1433116800
[919] 1435708800 1438387200 1441065600 1443657600 1446336000 1448928000 1451606400 1454284800 1456790400
[928] 1459468800
attr(,"tzone")
[1] "UTC"
attr(,"tclass")
[1] "Date"
$.indexCLASS
[1] "Date"
$tclass
[1] "Date"
$.indexTZ
[1] "UTC"
$tzone
[1] "UTC"
$src
[1] "FRED"
$updated
[1] "2016-05-18 21:30:19 JST"
$dimnames
$dimnames[[1]]
NULL
$dimnames[[2]]
[1] "PAYEMS"
$class
[1] "xts" "zoo"

> 

check attributes in XTS class

use attr function in order to check attributes in XTS class.

In the case below, all XTS timezone attribute is set to "UTC", while system timezone is NULL.

> attr(GDP_XTS,'tz')
[1] "UTC"
> attr(GDPC96,'tz')
[1] "UTC"
> attr(GDP,'tz')
[1] "UTC"
> attr(RRSFS,'tz')
[1] "UTC"
> Sys.getenv("TZ")
[1] ""

how to know current working directory

type getwd() on console.

> getwd()
[1] "/Users/<directory_name>/R_proj/tmp1"
> 

Please be aware that just type "getwd" doesn't work.

> getwd
function () 
.Internal(getwd())
<bytecode: 0x101fb2000>
<environment: namespace:base>

Please note that ,without any additional location information, the functions which read/write files look for the object in this directory. "read.csv" and "write.csv" are among them.

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

how to use mondate(1)

mondate　ライブラリを準備する。

>library(mondate)

関数を準備する。mycalという名前にする。名前衝突に注意すること。

mycal <- function() {
    y2 <- readline("year(yyyy)? ")
    m  <- readline("month(mm)? ")
    g <- trunc(MonthsBetween(as.mondate(as.Date("2016-01-01")),as.mondate(as.Date(ISOdate(y2,m,01)))))
    print(g)

}

how to add an entiry into XTS class.

> class(GDP_XTS)
[1] "xts" "zoo"
> last(GDP_XTS,n=4)[,1]
2015-10-01 16470.57
2015-11-01 16477.96
2015-12-01 16485.35
2016-01-01 16492.75

add an entry with index. this sample uses the index as "2016-02-01". data might be just garbage.

> last(append(GDP_XTS,c(as.Date("2016-02-01"))))[,1]
2016-02-01 2.145571e-314
> w <- append(GDP_XTS,c(as.Date("2016-02-01")))
> last(w)[,1]
2016-02-01 2.144668e-314

access entry by specifying index.

> w["2016-02-01"][,1]
2016-02-01 2.144668e-314
> w["2016-02-01"][,1] <- 16525

unable to insert a data entry where the index doesn't exit. Create the index before insert!

> w["2016-03-01"][,1] <- 16558
> last(w,n=4)[,1]
2015-11-01 16477.96
2015-12-01 16485.35
2016-01-01 16492.75
2016-02-01 16525.00
>