R TS VAR - Revision history

Rb: /* Model Estimation */

2016-09-16T10:09:38Z

‎Model Estimation

Rb: /* Introduction */

2016-09-16T10:08:24Z

‎Introduction

Rb: /* Introduction */

2016-09-16T10:08:02Z

‎Introduction

Rb at 23:44, 14 September 2016

2016-09-14T23:44:08Z

Rb: /* VAR Modeling */

2016-09-14T23:43:45Z

‎VAR Modeling

Rb: /* Model Estimation */

2016-09-14T23:40:25Z

‎Model Estimation

Rb: /* Impulse Response Functions (IRF) and re-ordering */

2016-09-14T23:37:23Z

‎Impulse Response Functions (IRF) and re-ordering

Rb: /* Transforming into (m)ts format */

2016-09-14T23:35:28Z

‎Transforming into (m)ts format

Rb: /* Some more resources */

2016-09-14T23:32:01Z

‎Some more resources

Rb: Created page with "= VAR Modeling = Ralf Becker = Introduction = Here we will demonstrate how to use R to model multivariate time series in a Vector Autoregressive Model (VAR). There are a v..."

2016-09-14T23:31:00Z

Created page with "= VAR Modeling = Ralf Becker = Introduction = Here we will demonstrate how to use R to model multivariate time series in a Vector Autoregressive Model (VAR). There are a v..."

New page

= VAR Modeling =

Ralf Becker

= Introduction =

Here we will demonstrate how to use R to model multivariate time series in a Vector Autoregressive Model (VAR).

There are a variety of ways to deal with time-series datasets in R. As discussed in another section of the [http://eclr.humanities.manchester.ac.uk/index.php/R_TimeSeries ECLR webpage] we could use the standard <code>ts</code> format.

In addition to treating the data as time-series objects we will need a package that does the hard lifting in terms of estimating the VAR, testing the model specification, producing forecasts or creating impulse response functions. Should these terms not make any sense to you then you should stop reading here and proceed by learning about VARs. Here we only show how to implement them in R. The package we use for this is the '''''vars''''' package (written by Bernard Pfaff) and hence, before you do anything else you shouldcheck whether you already have this package or whether you need to install it.

<pre class="r">install.packages("vars")</pre>
We then proceed by loading a few packages we need:

<pre class="r">setwd("T:/ECLR/R/VAR")
library(ggplot2)
library(xts)
library(reshape2) # this is required for the melt function below
library(vars)</pre>
= Import Data =

We use 10year maturity yields for the US, UK, German and French government bonds. The data come from the FRED database: [https://fred.stlouisfed.org/series/IRLTLT01GBM156N UK], [https://fred.stlouisfed.org/series/IRLTLT01USM156N US], [https://fred.stlouisfed.org/series/IRLTLT01DEM156N Germany] and [https://fred.stlouisfed.org/series/IRLTLT01FRM156N France].

We will use thge montly data from January 1975 to March 2015. These data are saved in a csv file with the four columns of data (LongRates.csv).

<pre class="r">rates <- read.csv("LongRates.csv")</pre>
Let's have a look at what these data look like:

<pre class="r">head(rates)</pre>
<pre>## DATE US Germany UK France
## 1 01/01/1975 7.50 9.4 15.07 10.90
## 2 01/02/1975 7.39 8.9 13.33 10.69
## 3 01/03/1975 7.73 8.9 12.24 10.34
## 4 01/04/1975 8.23 8.8 12.42 10.36
## 5 01/05/1975 8.06 8.6 13.11 10.31
## 6 01/06/1975 7.86 8.4 12.99 10.21</pre>
== Transforming into (m)ts format ==

Before we get cracking we want to ensure that the data are recognised as time-series data. Here we use the <code>ts</code> class which has also been used on the page for [http://eclr.humanities.manchester.ac.uk/index.php/R_TimeSeries univariate modelling].

We first pick each individual country's rate

<pre class="r">ger <- ts(rates$Germany, start=c(1975, 1), end=c(2015, 3), frequency=12)
fra <- ts(rates$France, start=c(1975, 1), end=c(2015, 3), frequency=12)
uk <- ts(rates$UK, start=c(1975, 1), end=c(2015, 3), frequency=12)
us <- ts(rates$US, start=c(1975, 1), end=c(2015, 3), frequency=12)
rates.ts = cbind(us,ger,uk,fra) </pre>
And let's look at a plot of these series:

<pre class="r">plot(rates.ts)</pre>
[[File:VARmodeling_files/figure-html/unnamed-chunk-6-1.png]]

In addition to the levels we shall also prepare an object with the differenced series.

<pre class="r">drates.ts <- diff(rates.ts)</pre>
And let's look at a plot of these series:

<pre class="r">plot(drates.ts)</pre>
[[File:VARmodeling_files/figure-html/unnamed-chunk-8-1.png]]

= Some preliminary data analysis =

Before you embark on a VAR analysis you may want to establish that the order of integration of your time series. You could, e.g., perform a unit root test for the German long rates. This can be done with the <code>ur.df</code> function which performs an Augmented Dickey-Fuller (ADF) test.

<pre class="r">summary(ur.df(rates.ts[, "ger"], type = "trend", lags = 2))</pre>
<pre>##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.96955 -0.12064 -0.00951 0.11437 0.86609
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1609140 0.0758616 2.121 0.0344 *
## z.lag.1 -0.0177337 0.0079895 -2.220 0.0269 *
## tt -0.0002992 0.0001306 -2.291 0.0224 *
## z.diff.lag1 0.3776531 0.0453528 8.327 8.92e-16 ***
## z.diff.lag2 -0.1055146 0.0454611 -2.321 0.0207 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1869 on 475 degrees of freedom
## Multiple R-squared: 0.1337, Adjusted R-squared: 0.1264
## F-statistic: 18.33 on 4 and 475 DF, p-value: 5.115e-14
##
##
## Value of test-statistic is: -2.2196 2.6479 2.7133
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36</pre>
If you know your unit root tests you can see that in this case we cannot reject the Null hypothesis of non-stationarity of this time series. You will find similar results for the other series.

Opinions on whether you should apply the VAR model to stationary or non-stationary data diverge somewhat with proponents on both sides. This is not the place for this debate. Here we decide to estimate a VAR model on stationary series and hence we will from now on use the differenced dataset <code>drates.ts</code>. (Of course you should confirm that the first differences are indeed stationary!)

= Estimating a VAR =

Now we demonstrate how to estimate the VAR in R.

== Lag selection ==

One of the first things you ought to check when estimating a VAR is how many lags you should use. We usually use the lag length that minimises some selection of information criterium. This is very easily done with the <code>VARselect</code> function. The options <code>lag.max = 8</code> and <code>type = "const"</code> indicate that we are looking for up to 8 lags and that we chose to include a constant as a deterministic variable into the VAR model. As usual use <code>?VARselect</code> to get details on the optional inputs.

<pre class="r">VARselect(drates.ts, lag.max = 8, type = "const")</pre>
<pre>## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 7 3 1 7
##
## $criteria
## 1 2 3 4
## AIC(n) -1.165290e+01 -1.172625e+01 -1.178755e+01 -1.177210e+01
## HQ(n) -1.158385e+01 -1.160196e+01 -1.160802e+01 -1.153732e+01
## SC(n) -1.147732e+01 -1.141021e+01 -1.133105e+01 -1.117513e+01
## FPE(n) 8.693850e-06 8.079051e-06 7.598975e-06 7.717867e-06
## 5 6 7 8
## AIC(n) -1.180622e+01 -1.179959e+01 -1.182123e+01 -1.181627e+01
## HQ(n) -1.151620e+01 -1.145433e+01 -1.142073e+01 -1.136053e+01
## SC(n) -1.106879e+01 -1.092170e+01 -1.080288e+01 -1.065746e+01
## FPE(n) 7.459759e-06 7.510543e-06 7.351413e-06 7.390077e-06</pre>
We can see that the different criteria select different lag lengths, e.g. the Hannan-Quinn Criterion (HQ) selects a lag length of 3, the Scwartz criterion (SC) selects a lag length of 1. Without any further consideration we shall from now on use lag length <code>p = 1</code>, mainly to reduce the size of the output in this illustration.

== Series ordering ==

For those of you in the know, you know that when using VARs (in particular when producing impulse response functions (see below) the order of the time-series plays an important role. Looking at the object <code>drates.ts</code> you will recognise that the order is <code>us, ger, uk, fra</code>.

We will later see how to re-order these.

== Model Estimation ==

The function that does the heavy liftung is the <code>VAR</code> function. Let's estimate a VAR of order 3 (<code>p = 3</code>) for the rate differences, including a constant (<code>type = "const"</code>), but no trend into the VAR.

<pre class="r">mod1 <- VAR(drates.ts, p = 1, type = "const")</pre>
This saves the estimated model in the object <code>mod1</code>, and we can show some basic output using the <code>summary(mod1)</code> command. The output below shows all estimated parameters in all 4 equations as well as covariance and correlation matrices of the estimated residuals.

== Specification testing ==

We will often be concerned whether the chosen lag length was sufficiently large to capture all the systematic temporal correlation in the time-series. In case it is not we should expect to find a certain degree of autocorrelation in the estimated residuals.

So, we need to find to know how to test for resdual autocorrelation. Please use the help function to read about the parameter choices. <code>lags.pt = 12</code> selects the maximum order of autocorrelation that is being tested for.

<pre class="r">ser12 <- serial.test(mod1, lags.pt = 12, type = "PT.asymptotic")
print(ser12)</pre>
<pre>##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object mod1
## Chi-squared = 380.75, df = 176, p-value < 2.2e-16</pre>
<pre>## $serial
##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object mod1
## Chi-squared = 380.75, df = 176, p-value < 2.2e-16</pre>
The p-value is very small indicating that the null hypothesis of no autocorrelation is rejected. This is a signal that we really ought to increase tha lag length. Incidentally, it turns out that only once you choose VAR orders of 7 or higher residual autocorrelation is by and large eliminated. Despite this we shall continue working with the very basic VAR(1) model.

== Impulse Response Functions (IRF) and re-ordering ==

One way in which VAR models are used is by interpreting the impulse response functions. To do this we use the <code>irf</code> function. When looking at the output of the estimated model <code>mod1</code> we could see that the estimated residuals are quite strongly correlated with each other. This means that it is very important to generate orthogonal impulse response functions. The standard way to do that is by employing a Cholesky decomposition (option <code>ortho = TRUE</code> in the <code>irf</code> function).

<pre class="r">irf1 <- irf(mod1, impulse = "ger", response = c("us", "ger", "uk", "fra"), ortho = TRUE, n.ahead = 6)</pre>
In this function call the option <code>impulse = "ger"</code> specifies that we are looking at German shocks and <code>response = c("us", "ger", "uk", "fra")</code> indicates for which series we want to see the IRF.

This saves the impulse response function in <code>irf1</code>. You can either print the values by <code>print(irf1)</code>, or, more pleasingly plot them.

<pre class="r">plot(irf1)</pre>
[[File:VARmodeling_files/figure-html/unnamed-chunk-14-1.png]]

Here we can see that the strongest impact of a German shock is to the German long rate. UK and France rates also react but not as strongly. The US rate doesn't reacht immeadiately and barely after that. The fact that the US rate doesn't react to a German shock in the period of the shock (period 0) is by design. Recall that the ordering of the series is (US, GER, UK, FRA) implying that in the shock period a German shock cannot impact on the US rate.

Of course we may want to change the ordering and re-evaluate the IRF. In the next block of code we first re-order (change the position of US and GER), then re-estimate and present the IRF for the new ordering.

<pre class="r">drates.ts <- drates.ts[,c("ger","us","uk","fra")]
mod1r <- VAR(drates.ts, p = 1, type = "const")
irf1r <- irf(mod1, impulse = "ger", response = c("us", "ger", "uk", "fra"), ortho = TRUE, n.ahead = 6)
plot(irf1r)</pre>
[[File:VARmodeling_files/figure-html/unnamed-chunk-15-1.png]]

As it turns out, even though now the US could react instantaneously to a GER shock, it doesn't. This indicates that the ordering of th US at the top of the heap is justified in this instance.

= Some more resources =

A few good graphing resources for ggplot are here:

* A detailed document on the <code>vars</code> package by its creator Bernard Pfaff can be found [https://cran.r-project.org/web/packages/vars/vignettes/vars.pdf here].

@@ Line 135: / Line 135: @@
 == Model Estimation ==
-The function that does the heavy liftung is the <code>VAR</code> function. Let's estimate a VAR of order 3 (<code>p = 3</code>) for the rate differences, including a constant (<code>type = &quot;const&quot;</code>), but no trend into the VAR.
+The function that does the heavy lifting is the <code>VAR</code> function. Let's estimate a VAR of order 3 (<code>p = 3</code>) for the rate differences, including a constant (<code>type = &quot;const&quot;</code>), but no trend into the VAR.
 <pre class="r">mod1 &lt;- VAR(drates.ts, p = 1, type = &quot;const&quot;)</pre>

← Older revision		Revision as of 23:44, 14 September 2016
Line 1:		Line 1:
−	~~= VAR Modeling =~~
−	~~This page introduces the ways in which to implement a VAR model in R.~~
−
	= Introduction =		= Introduction =

@@ Line 10: / Line 10: @@
 We then proceed by loading a few packages we need:
-<pre class="r">setwd(&quot;T:/ECLR/R/VAR&quot;)
+<pre class="r">setwd(&quot;T:/ECLR/R/VAR&quot;) # replace with your working directory
 library(xts)
 library(vars)</pre>

@@ Line 1: / Line 1: @@
 = VAR Modeling =
+This page introduces the ways in which to implement a VAR model in R.
 = Introduction =

@@ Line 144: / Line 144: @@
 <pre class="r">mod1 &lt;- VAR(drates.ts, p = 1, type = &quot;const&quot;)</pre>
 This saves the estimated model in the object <code>mod1</code>, and we can show some basic output using the <code>summary(mod1)</code> command. The output below shows all estimated parameters in all 4 equations as well as covariance and correlation matrices of the estimated residuals.
 == Specification testing ==

@@ Line 176: / Line 176: @@
 <pre class="r">plot(irf1)</pre>
-[[File:VARmodeling_files/figure-html/unnamed-chunk-14-1.png]]<!-- -->
+[[File:IRF1.png]]<!-- -->
 Here we can see that the strongest impact of a German shock is to the German long rate. UK and France rates also react but not as strongly. The US rate doesn't reacht immeadiately and barely after that. The fact that the US rate doesn't react to a German shock in the period of the shock (period 0) is by design. Recall that the ordering of the series is (US, GER, UK, FRA) implying that in the shock period a German shock cannot impact on the US rate.
@@ Line 186: / Line 186: @@
 irf1r &lt;- irf(mod1, impulse = &quot;ger&quot;, response = c(&quot;us&quot;, &quot;ger&quot;, &quot;uk&quot;, &quot;fra&quot;), ortho = TRUE, n.ahead = 6)
 plot(irf1r)</pre>
-[[File:VARmodeling_files/figure-html/unnamed-chunk-15-1.png]]<!-- -->
+[[File:IRF2.png]]<!-- -->
 As it turns out, even though now the US could react instantaneously to a GER shock, it doesn't. This indicates that the ordering of th US at the top of the heap is justified in this instance.

@@ Line 50: / Line 50: @@
 <pre class="r">plot(rates.ts)</pre>
-[[File:VARmodeling_files/figure-html/unnamed-chunk-6-1.png]]<!-- -->
+[[File:rates.png]]<!-- -->
 In addition to the levels we shall also prepare an object with the differenced series.
@@ Line 58: / Line 58: @@
 <pre class="r">plot(drates.ts)</pre>
-[[File:VARmodeling_files/figure-html/unnamed-chunk-8-1.png]]<!-- -->
+[[File:drates.png]]<!-- -->
 = Some preliminary data analysis =

@@ Line 192: / Line 192: @@
 = Some more resources =
-A few good graphing resources for ggplot are here:
+A few good resources are here:
 * A detailed document on the <code>vars</code> package by its creator Bernard Pfaff can be found [https://cran.r-project.org/web/packages/vars/vignettes/vars.pdf here].