Analysis of Integrated and Cointegrated Time Series with R (Pfaff 2008 Springer)
There are two main sections: brief explanatory text, and code provided with the book.
1 Chapter intro and summary
The book can be found at: https://sci-hub.tw/10.1007/978-0-387-75967-8
Here the text is quoted literally without edits. For each chapter the introductury paragraph and the summary are shown, to provide context for the code.
1.1 Univariate
1.1.1 Intro
Although this book has integration and cointegration analysis as its theme, it is nevertheless a necessity to first introduce some concepts of stochastic processes as well as the stationary ARMA model class. Having paved this route, the next steps (i.e., the introduction of nonstationary, unit root, and long-memory processes) will follow in Chapter 3.
1.1.2 Summary
In this first chapter, the definition of a time series and the concept of its datagenerating process have been introduced. You should now be familiar with how to characterize a time series by its moments and distinguish the different concepts of stationarity. Two model classes for a time series have been introduced, namely the autoregressive and the moving average models, as well as a combination thereof. You should be able to detect and distinguish the order of these models by investigating its autocorrelation and partial auto-correlation functions. Finally, the Box-Jenkins approach to time series analysis has been presented. It is decomposed into three stages: specification of a tentative model order, estimation, and diagnostic checking. So far, we have restricted the presentation to stationary time series only. At first sight, this focus might seem to be too myopic given that many time series cannot be characterized by a stationary process, in particular in macroeconomic and financial data sets. Therefore, in Chapter 3, non-stationary time series processes and how they can be transformed to achieve stationarity are discussed
1.2 Multivariate
1.2.1 Intro
This is the second chapter that presents models confined to stationary time series, but now in the context of multivariate analysis. Vector autoregressive models and structural vector autoregressive models are introduced. The analytical tools of impulse response functions, forecast error variance decomposition, and Granger causality, as well as forecasting and diagnostic tests, are outlined. As will be shown later, these concepts can be applied to cointegrated systems, too.
1.2.2 Summary
In this chapter, the analysis of stationary time series has been extended to multivariate models and their associated statistical tests and methods. In particular, VAR- and SVAR-models have been introduced, where the former can be interpreted as the reduced-form counterparts of SVAR-models. Both model classes have been illustrated by artificial data sets. It has been outlined how a suitable lag length can be empirically determined and what kind of diagnostic tests are at hand for checking the assumptions about the multivariate error process. The different concepts of causality analysis and forecasting with VAR-models have been shown. For investigating the dynamic interactions between variables, the impulse response functions and forecast error variance decomposition have been introduced. These tools are implemented as methods for VAR- and SVAR-models alike. The results can be obtained and plotted swiftly with the functions included in package vars. An overview of the package’s structure is presented in Table 2.9.
1.3 Nonstationary
1.3.1 Intro
In this chapter, models for non-stationary time series are introduced. Before the characteristics of unit processes are presented, the differences between trend- and difference-stationary models are outlined. In the last section, long-memory processes (i.e., fractionally integrated processes) are presented as a bridge between stationary and unit root processes.
1.3.2 Summary
In this chapter, a more encompassing data-generating process that was introduced into the literature by Campbell and Perron [1991] has been presented. You should now be familiar with the concepts of trend- versus differencestationary and the decomposition of a time series into a deterministic trend, a stochastic trend, and a cyclical component. Furthermore, unit root processes have been introduced as a subclass of random walk processes. How one applies a sequential testing strategy to detect the underlying data-generating process of a possible non-stationary time series was discussed. The important definitions of integrated, seasonally integrated, and fractionally integrated time series processes have been presented, too, where the latter can be viewed as a bridge between stationary and unit root processes, thereby closing the circle of the exposition in the first two chapters. So far, we have addressed univariate time series analysis and multivariate analysis in the context of stationary VAR and SVAR models only. The obstacles and solutions for multivariate models with non-stationary data are the subject of the next and last chapter of Part I.
1.4 Cointegration
1.4.1 Intro
In the previous chapters, a brief explanation of univariate and multivariate time series models and their characteristics was presented. The focus of this chapter is on the simultaneous modeling of time series and inferences of the relationships between them if some or all of them are integrated processes of order one. As will be shown, the degree of integration and a careful examination of the data-generating processes are of utmost importance. We will begin by briefly reviewing the case of a spurious regression before we proceed by providing a definition of cointegration and its error-correction representation. In the last section, the more encompassing vector error-correction models are presented.
1.4.2 Summary
By outlining the spurious regression problem in the first section, you should have been alerted to the pitfalls when integrated time series are modeled in a multivariate context; recall the rule of thumb for detecting such nonsense regressions. In Section 4.2, the solution to this problem was presented, namely the definition of cointegration as a linear combination with a degree of lower integratedness than the two integrated processes to be investigated. In this respect, it was pointed out first that if two time series are cointegrated, then an error-correction mechanism exists and vice versa and, second, that in the case of two cointegrated I(1)-variables, Granger causality must exist in at least one direction. An empirical example of these issues is presented in Section 7.1. The shortcoming of the Engle-Granger procedure is that in the case of more than two integrated variables not only one cointegrating vector can exist. In fact, by applying the Engle-Granger two-step procedure in cases with more than one cointegrating vector, one would estimate a linear combination of these vectors. An encompassing definition of cointegration and the model class of VECM has been introduced to cope adequately with such instances. It has been shown that two forms of a VECM exist and that the inference with regard to the order of the space spanned by the linearly independent cointegrating vectors is the same. You should recall that neither the cointegrating vectors nor the adjustment matrix can uniquely be determined. Instead a normalization of one element of β to one is applied. Up to now, the likelihood-based inference in cointegrated vector autoregressive models has been confined to determining the cointegration rank only. Testing various restrictions placed on the cointegration vectors and the loading matrix as well as a combination thereof are presented in Subsections 8.1.3 and 8.1.4.
1.5 Order of Integration
1.5.1 Intro
This chapter is the first in which the theoretical aspects laid out in Part I of the book are put into “practice.” We begin by introducing the most commonly employed unit root tests in econometrics: the Dickey-Fuller test and its extensions. To discriminate between trend- and difference-stationary time series processes, a sequential testing strategy is described. Other unit root tests encountered in applied research are presented in the ensuing sections.
1.5.2 Summary
In this first chapter of Part II, various unit root tests have been applied to real data sets. The sequential testing strategy of the ADF test outlined in Section 3.2 has been applied to U.K. consumption. Because the data-generating process is unknown, it is recommended to go through these steps rather than merely apply one-test regressions as in Equations (5.1a)–(5.1c). Furthermore, a spherical error term should always be ensured by supplying sufficient lagged endogenous variables. Next, the Phillips-Perron test has been applied to the same data set. In principle, the difference between the two tests is that the latter uses a non-parametric correction that captures weak dependence and heterogeneity of the error process. As pointed out in Section 5.3, the relatively low power of both tests due to the fact that a unit root process is specified as the null hypothesis must be considered as a shortcoming. The ERS tests ameliorate this problem and should therefore be preferred. Furthermore, the nuisance parameters have different interpretations if either the null or the alternative hypothesis is true. The SP test addresses this problem explicitly and allows inclusion of higher polynomials in the deterministic part. However, all tests suffer from an ill-specified null hypothesis. The KPSS test, as a test for stationarity, correctly addresses the hypothesis specification from the viewpoint of conservative testing. Anyway, unfortunately there is no clear-cut answer to the question of which test should be applied to a data set. A combination of some of the above-mentioned tests with the inclusion of opposing null hypotheses therefore seems to be a pragmatic approach in practice
1.6 Further Tests
1.6.1 Intro
In Chapter 5, various unit root tests were introduced and compared with each other. This chapter deals with two further topics. First, the case of structural breaks in a time series and how this affects the inference about the degree of integratedness is considered. Second, the issue of seasonal unit roots is discussed, as it was only briefly touched on in Section 3.2.
1.6.2 Summary
In this chapter, the analysis of integrated time series has been amended in two important ways. First, it has been shown that in the case of a structural break the test conclusion about the presence of a unit root in a time series can be biased toward accepting it. Therefore, if a priori knowledge of a structural shift exists or a break in the series is evident by eye-spotting, one should either use the Perron or the Zivot and Andrews test, respectively. Second, if the correlogram gives hindsight of seasonality in the time series, one should apply a seasonal unit root test. A complete analysis of a possibly integrated time series would therefore begin by testing whether breaks and/or stochastic seasonality exist, and depending on this outcome, unit root tests should be applied as shown in Chapter 5. After all, when the null hypothesis of a unit root must be rejected, it should be checked whether long-memory behavior is present as shown in Section 3.3.
1.7 Single-Equation
1.7.1 Intro
This is the first chapter of the third and last part of this book. The cointegration methodology is first presented for the case of single-equation models. The Engle-Granger two-step procedure is demonstrated by estimating a consumption function and its error-correction form for the United Kingdom as in Holden and Perman [1994]. In the Section 7.2, the method proposed by Phillips and Ouliaris [1990] is applied to the same data set. The application and inferences of a vector error-correction model are saved for Chapter 8.
1.7.2 Summary
In this first chapter of Part III, two single-equation methods have been presented. The advantage of the Engle-Granger two-step procedure is its ease of implementation. However, the results are dependent on how the long-run equation is specified. In most cases, it might be obvious which variable enters on the left-hand side of the equation; i.e., to which variable the cointegrating vector should be normalized. Unfortunately, this is only true in most cases, and, as anecdotal evidence, an income function rather than a consumption function could have been specified as an ECM in R code 7.2. It is therefore advisable to employ the cointegration test of Phillips and Ouliaris, which is irrelevant to normalization. As mentioned, the insights gained with respect to the cointegrating relationship are limited in the case of more than two variables. The next chapter is therefore dedicated to the inference in cointegrated systems.
1.8 Multiple-Equation
1.8.1 Intro
In this chapter, the powerful tool of likelihood-based inference in cointegrated vector autoregressive models (VECMs) is discussed. In the first section, the specification and assumptions of a VECM are introduced. In the following sections, the problems of determining the cointegration rank, testing for weak exogenity, and testing of various restrictions placed on the cointegrating vectors are discussed. The topic of VECMs that are contaminated by a one-time structural shift and how this kind of model can be estimated are presented. This chapter concludes with an exposition of structural vector error correction models.
1.8.2 Summary
In this last chapter of the book, likelihood-based inference in cointegrated vector autoregressive models has been presented. It has been shown how to determine the cointegration rank and, depending on that outcome, how to specify and test the validity of restrictions placed on the cointegrating and the weighting matrices. This methodology offers the researcher a powerful tool to investigate the relationships in a system of cointegrated variables more thoroughly compared with the single-equation methods presented in Chapter 7. Furthermore, it has been shown how one can employ this methodology in light of a structural shift at an unknown point in time. The chapter concluded with an exposition of the SVEC analysis applied to a macroeconomic data set for Canada. ## Appendix
2 Code
Source is at https://www.pfaffikus.de/books/spex2/
2.1 Univariate
2.1.1 Simulation of AR(1)-process with φ = 0.9
#title: "Chunk 1-1"
set.seed(123456)
y <- stats::arima.sim(n = 100, list(ar = 0.9), innov = rnorm(100))
op <- graphics::par(no.readonly = TRUE)
graphics::layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE))
stats::plot.ts(y, ylab = '')
stats::acf(
y,
main = 'Autocorrelations',
ylab = '',
ylim = c(-1, 1),
ci.col = "black"
)
stats::pacf(
y,
main = 'Partial Autocorrelations',
ylab = '',
ylim = c(-1, 1),
ci.col = "black"
)
graphics::par(op) #-fig1.1
# library(fArma)
# armaSim()
# [fArma is a discontinued package]2.1.2 Estimation of AR(2)-process with φ1 = 0.6 and φ2 = −0.28
series <- rnorm(1000)
y.st <- filter(series, filter = c(0.6,-0.28),
method = 'recursive')
ar2.st <- arima(
y.st,
c(2, 0, 0),
include.mean = FALSE,
transform.pars = FALSE,
method = "ML"
)
ar2.st$coef## ar1 ar2
## 0.6147474 -0.2835138polyroot(c(1,-ar2.st$coef))## [1] 1.084158+1.533547i 1.084158-1.533547iMod(polyroot(c(1,-ar2.st$coef)))## [1] 1.878075 1.878075root.comp <- Im(polyroot(c(1,-ar2.st$coef)))
root.real <- Re(polyroot(c(1,-ar2.st$coef)))
# Plotting the roots in a unit circle
x <- seq(-1, 1, length = 1000)
y1 <- sqrt(1 - x ^ 2)
y2 <- -sqrt(1 - x ^ 2)
plot( #fig 1.5
c(x, x),
c(y1, y2),
xlab = 'Real part',
ylab = 'Complex part',
type = 'l',
main = 'Unit Circle',
ylim = c(-2, 2),
xlim = c(-2, 2)
)
abline(h = 0)
abline(v = 0)
points(Re(polyroot(c(1,-ar2.st$coef))),
Im(polyroot(c(1,-ar2.st$coef))), pch = 19)
legend(-1.5,-1.5, legend = "Roots of AR(2)", pch = 19)
#-----comment
#fig 1.4 is not done so copy code from 101
y <- y.st
op <- graphics::par(no.readonly = TRUE)
graphics::layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE))
stats::plot.ts(y, ylab = '')
stats::acf(
y,
main = 'Autocorrelations',
ylab = '',
ylim = c(-1, 1),
ci.col = "black"
)
stats::pacf(
y,
main = 'Partial Autocorrelations',
ylab = '',
ylim = c(-1, 1),
ci.col = "black"
)
graphics::par(op) #-fig1.42.1.3 Box-Jenkins: U.S. unemployment rate
library(urca)
data(npext)
y <- ts(
na.omit(npext$unemploy),
start = 1909,
end = 1988,
frequency = 1
)
op <- par(no.readonly = TRUE) #fig 1.6
layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE))
plot(y, ylab = "unemployment rate (logarithm)")
acf(y,
main = 'Autocorrelations',
ylab = '',
ylim = c(-1, 1))
pacf(y,
main = 'Partial Autocorrelations',
ylab = '',
ylim = c(-1, 1))
par(op)
## tentative ARMA(2,0)
arma20 <- arima(y, order = c(2, 0, 0))
ll20 <- logLik(arma20)
aic20 <- arma20$aic
res20 <- residuals(arma20)
Box.test(res20, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: res20
## X-squared = 21.914, df = 20, p-value = 0.3452shapiro.test(res20)##
## Shapiro-Wilk normality test
##
## data: res20
## W = 0.99313, p-value = 0.9501## alternative specifications
## ARMA(3,0)
arma30 <- arima(y, order = c(3, 0, 0))
ll30 <- logLik(arma30)
aic30 <- arma30$aic
lrtest <- as.numeric(2 * (ll30 - ll20))
chi.pval <- pchisq(lrtest, df = 1, lower.tail = FALSE)
## ARMA(1,1)
arma11 <- arima(y, order = c(1, 0, 1))
ll11 <- logLik(arma11)
aic11 <- arma11$aic
tsdiag(arma11) #fig 1.8 
res11 <- residuals(arma11)
Box.test(res11, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: res11
## X-squared = 15.14, df = 20, p-value = 0.7683shapiro.test(res11)##
## Shapiro-Wilk normality test
##
## data: res11
## W = 0.98617, p-value = 0.5456## Using auto.arima()
library(forecast)## Warning: package 'forecast' was built under R version 3.6.3## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooauto.arima(
y,
max.p = 3,
max.q = 3,
start.p = 1,
start.q = 1,
ic = "aic"
)## Series: y
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.5272 0.5487 1.6934
## s.e. 0.1221 0.1456 0.1546
##
## sigma^2 estimated as 0.1917: log likelihood=-46.51
## AIC=101.01 AICc=101.55 BIC=110.542.1.4 Box-Jenkins: Predictions of the U.S. unemployment rate
## Forecasts
arma11.pred <- predict(arma11, n.ahead = 10)
predict <- ts(c(rep(NA, length(y) - 1), y[length(y)],
arma11.pred$pred),
start = 1909,
frequency = 1)
upper <- ts(
c(rep(NA, length(y) - 1), y[length(y)],
arma11.pred$pred + 2 * arma11.pred$se),
start = 1909,
frequency = 1
)
lower <- ts(
c(rep(NA, length(y) - 1), y[length(y)],
arma11.pred$pred - 2 * arma11.pred$se),
start = 1909,
frequency = 1
)
observed <- ts(c(y, rep(NA, 10)), start = 1909,
frequency = 1)
## Plot of actual and forecasted values
plot(observed, # fig 1.9
type = "l",
ylab = "Actual and predicted values",
xlab = "")
lines(predict, col = "blue", lty = 2)
lines(lower, col = "red", lty = 5)
lines(upper, col = "red", lty = 5)
abline(v = 1988, col = "gray", lty = 3)
2.2 Multivariate
2.2.1 Simulation of VAR(2)-process
## Simulate VAR(2)-data
library(dse)## Warning: package 'dse' was built under R version 3.6.3## Loading required package: tfplot## Warning: package 'tfplot' was built under R version 3.6.3## Loading required package: tframe##
## Attaching package: 'dse'## The following objects are masked from 'package:forecast':
##
## forecast, is.forecast## The following objects are masked from 'package:stats':
##
## acf, simulatelibrary(vars)## Loading required package: MASS## Loading required package: strucchange## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Loading required package: lmtest##
## Attaching package: 'vars'## The following objects are masked from 'package:dse':
##
## roots, stability## Setting the lag-polynomial A(L)
Apoly <- array(c(1.0,-0.5, 0.3, 0,
0.2, 0.1, 0,-0.2,
0.7, 1, 0.5,-0.3) ,
c(3, 2, 2))
## Setting Covariance to identity-matrix
B <- diag(2)
## Setting constant term to 5 and 10
TRD <- c(5, 10)
## Generating the VAR(2) model
var2 <- dse::ARMA(A = Apoly, B = B, TREND = TRD)
## Simulating 500 observations
varsim <- dse::simulate(
var2,
sampleT = 500,
noise = list(w = matrix(
rnorm(1000),
nrow = 500, ncol = 2
)),
rng = list(seed = c(123456))
)
## Obtaining the generated series
vardat <- matrix(varsim$output, nrow = 500, ncol = 2)
colnames(vardat) <- c("y1", "y2")
## Plotting the series
plot.ts(vardat, main = "", xlab = "") #fig 2.1
## Determining an appropriate lag-order
infocrit <- VARselect(vardat, lag.max = 3,
type = "const")
## Estimating the model
varsimest <- VAR(
vardat,
p = 2,
type = "const",
season = NULL,
exogen = NULL
)
## Alternatively, selection according to AIC
varsimest <- VAR(vardat,
type = "const",
lag.max = 3,
ic = "SC")
## Checking the roots
roots <- roots(varsimest)
#-add
print(varsimest) # for tables 2.1 & 2.2, approx correct##
## VAR Estimation Results:
## =======================
##
## Estimated coefficients for equation y1:
## =======================================
## Call:
## y1 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## y1.l1 y2.l1 y1.l2 y2.l2 const
## 0.5188174 0.1675959 -0.3319501 -0.7591715 5.7608039
##
##
## Estimated coefficients for equation y2:
## =======================================
## Call:
## y2 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## y1.l1 y2.l1 y1.l2 y2.l2 const
## -0.2026751 -0.3985342 -0.1102268 0.3238735 8.9659360print(roots) #table 2.3 'eigenvalues', approx correct## [1] 0.7757285 0.6588502 0.6116253 0.61162532.2.2 Diagnostic tests of VAR(2)-process
## testing serial correlation
args(serial.test)## function (x, lags.pt = 16, lags.bg = 5, type = c("PT.asymptotic",
## "PT.adjusted", "BG", "ES"))
## NULL## Portmanteau-Test
var2c.serial <- serial.test(varsimest, lags.pt = 16,
type = "PT.asymptotic")
var2c.serial##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 49.741, df = 56, p-value = 0.7093plot(var2c.serial, names = "y1") #fig 2.2 [not identical but similar]
plot(var2c.serial, names = "y2") #fig 2.3 [not identical but similar]
## testing heteroscedasticity
args(arch.test)## function (x, lags.single = 16, lags.multi = 5, multivariate.only = TRUE)
## NULLvar2c.arch <- arch.test(varsimest, lags.multi = 5,
multivariate.only = TRUE)
var2c.arch##
## ARCH (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 46.494, df = 45, p-value = 0.4106## testing for normality
args(normality.test)## function (x, multivariate.only = TRUE)
## NULLvar2c.norm <- normality.test(varsimest,
multivariate.only = TRUE)
var2c.norm## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 5.9826, df = 4, p-value = 0.2004
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 3.6933, df = 2, p-value = 0.1578
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 2.2894, df = 2, p-value = 0.3183## class and methods for diganostic tests
class(var2c.serial)## [1] "varcheck"class(var2c.arch)## [1] "varcheck"class(var2c.norm)## [1] "varcheck"methods(class = "varcheck")## [1] plot print
## see '?methods' for accessing help and source code## Plot of objects "varcheck"
args(vars:::plot.varcheck)## function (x, names = NULL, main.resid = NULL, main.hist = NULL,
## main.acf = NULL, main.pacf = NULL, main.acf2 = NULL, main.pacf2 = NULL,
## ylim.resid = NULL, ylim.hist = NULL, ylab.resid = NULL, xlab.resid = NULL,
## xlab.acf = NULL, lty.resid = NULL, lwd.resid = NULL, col.resid = NULL,
## col.edf = NULL, lag.acf = NULL, lag.pacf = NULL, lag.acf2 = NULL,
## lag.pacf2 = NULL, mar = par("mar"), oma = par("oma"), ...)
## NULLplot(var2c.serial, names = "y1")
#add
print(var2c.serial) #table 2.4 portmanteau##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 49.741, df = 56, p-value = 0.7093## $serial
##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 49.741, df = 56, p-value = 0.7093print(var2c.arch) #table 2.4 ARCH##
## ARCH (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 46.494, df = 45, p-value = 0.4106## $arch.mul
##
## ARCH (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 46.494, df = 45, p-value = 0.4106print(var2c.norm) #table 2.4 jarque-bera, kurtosis only, skew only## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 5.9826, df = 4, p-value = 0.2004
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 3.6933, df = 2, p-value = 0.1578
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 2.2894, df = 2, p-value = 0.3183## $jb.mul
## $jb.mul$JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 5.9826, df = 4, p-value = 0.2004
##
##
## $jb.mul$Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 3.6933, df = 2, p-value = 0.1578
##
##
## $jb.mul$Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object varsimest
## Chi-squared = 2.2894, df = 2, p-value = 0.31832.2.3 Empirical fluctuation processes
reccusum <- stability(varsimest,
type = "OLS-CUSUM")
fluctuation <- stability(varsimest,
type = "fluctuation")
#-add to get plots
plot(reccusum) #fig 2.4
plot(fluctuation) #fig 2.5
2.2.4 Causality analysis of VAR(2)-process
## Causality tests
## Granger and instantaneous causality
var.causal <- causality(varsimest, cause = "y2")
#add
print(var.causal) #table 2.5## $Granger
##
## Granger causality H0: y2 do not Granger-cause y1
##
## data: VAR object varsimest
## F-Test = 212.69, df1 = 2, df2 = 986, p-value < 2.2e-16
##
##
## $Instant
##
## H0: No instantaneous causality between: y2 and y1
##
## data: VAR object varsimest
## Chi-squared = 0.016052, df = 1, p-value = 0.89922.2.5 Forecasts of VAR-process
## Forecasting objects of class varest
args(vars:::predict.varest)## function (object, ..., n.ahead = 10, ci = 0.95, dumvar = NULL)
## NULLpredictions <- predict(varsimest, n.ahead = 25,
ci = 0.95)
class(predictions)## [1] "varprd"args(vars:::plot.varprd)## function (x, plot.type = c("multiple", "single"), names = NULL,
## main = NULL, col = NULL, lty = NULL, lwd = NULL, ylim = NULL,
## ylab = NULL, xlab = NULL, nc, mar = par("mar"), oma = par("oma"),
## ...)
## NULL## Plot of predictions for y1
plot(predictions, names = "y1") #fig 2.6
## Fanchart for y2
args(fanchart)## function (x, colors = NULL, cis = NULL, names = NULL, main = NULL,
## ylab = NULL, xlab = NULL, col.y = NULL, nc, plot.type = c("multiple",
## "single"), mar = par("mar"), oma = par("oma"), ...)
## NULLfanchart(predictions, names = "y2") #fig 2.7
2.2.6 IRA of VAR-process
## Impulse response analysis
irf.y1 <- irf(varsimest, impulse = "y1", #fig 2.8
response = "y2", n.ahead = 10,
ortho = FALSE, cumulative = FALSE,
boot = FALSE, seed = 12345)
args(vars:::plot.varirf)## function (x, plot.type = c("multiple", "single"), names = NULL,
## main = NULL, sub = NULL, lty = NULL, lwd = NULL, col = NULL,
## ylim = NULL, ylab = NULL, xlab = NULL, nc, mar.multi = c(0,
## 4, 0, 4), oma.multi = c(6, 4, 6, 4), adj.mtext = NA,
## padj.mtext = NA, col.mtext = NA, ...)
## NULLplot(irf.y1)
irf.y2 <- irf(varsimest, impulse = "y2", #fig 2.9
response = "y1", n.ahead = 10,
ortho = TRUE, cumulative = TRUE,
boot = FALSE, seed = 12345)
plot(irf.y2)
2.2.7 FEVD of VAR-process
## Forecast error variance decomposition
fevd.var2 <- fevd(varsimest, n.ahead = 10)
args(vars:::plot.varfevd)## function (x, plot.type = c("multiple", "single"), names = NULL,
## main = NULL, col = NULL, ylim = NULL, ylab = NULL, xlab = NULL,
## legend = NULL, names.arg = NULL, nc, mar = par("mar"), oma = par("oma"),
## addbars = 1, ...)
## NULLplot(fevd.var2, addbars = 2) #fig 2.10
2.2.8 SVAR: A-model
library(dse)
library(vars)
## A-model
Apoly <- array(c(1.0, -0.5, 0.3, 0.8,
0.2, 0.1, -0.7, -0.2,
0.7, 1, 0.5, -0.3) ,
c(3, 2, 2))
## Setting covariance to identity-matrix
B <- diag(2)
## Generating the VAR(2) model
svarA <- ARMA(A = Apoly, B = B)
## Simulating 500 observations
svarsim <- simulate(svarA, sampleT = 500,
rng = list(seed = c(123456)))
## Obtaining the generated series
svardat <- matrix(svarsim$output, nrow = 500, ncol = 2)
colnames(svardat) <- c("y1", "y2")
## Estimating the VAR
varest <- VAR(svardat, p = 2, type = "none")
## Setting up matrices for A-model
Amat <- diag(2)
Amat[2, 1] <- NA
Amat[1, 2] <- NA
## Estimating the SVAR A-type by direct maximisation
## of the log-likelihood
args(SVAR)## function (x, estmethod = c("scoring", "direct"), Amat = NULL,
## Bmat = NULL, start = NULL, max.iter = 100, conv.crit = 1e-07,
## maxls = 1, lrtest = TRUE, ...)
## NULLsvar.A <- SVAR(varest, estmethod = "direct",
Amat = Amat, hessian = TRUE)
#add
print(svar.A) #tale 2.6##
## SVAR Estimation Results:
## ========================
##
##
## Estimated A matrix:
## y1 y2
## y1 1.0000 -0.6975
## y2 0.8571 1.00002.2.9 SVAR: B-model
library(dse)
library(vars)
## B-model
Apoly <- array(c(1.0, -0.5, 0.3, 0,
0.2, 0.1, 0, -0.2,
0.7, 1, 0.5, -0.3) ,
c(3, 2, 2))
## Setting covariance to identity-matrix
B <- diag(2)
B[2, 1] <- -0.8
## Generating the VAR(2) model
svarB <- ARMA(A = Apoly, B = B)
## Simulating 500 observations
svarsim <- simulate(svarB, sampleT = 500,
rng = list(seed = c(123456)))
svardat <- matrix(svarsim$output, nrow = 500, ncol = 2)
colnames(svardat) <- c("y1", "y2")
varest <- VAR(svardat, p = 2, type = "none")
## Estimating the SVAR B-type by scoring algorithm
## Setting up the restriction matrix and vector
## for B-model
Bmat <- diag(2)
Bmat[2, 1] <- NA
svar.B <- SVAR(varest, estmethod = "scoring",
Bmat = Bmat, max.iter = 200)
#add
print(svar.B) #table 2.7##
## SVAR Estimation Results:
## ========================
##
##
## Estimated B matrix:
## y1 y2
## y1 1.0000 0
## y2 -0.8439 12.2.10 SVAR: Impulse response analysis
## Impulse response analysis of SVAR A-type model
args(vars:::irf.svarest)## function (x, impulse = NULL, response = NULL, n.ahead = 10, ortho = TRUE,
## cumulative = FALSE, boot = TRUE, ci = 0.95, runs = 100, seed = NULL,
## ...)
## NULLirf.svara <- irf(svar.A, impulse = "y1",
response = "y2", boot = FALSE)
args(vars:::plot.varirf)## function (x, plot.type = c("multiple", "single"), names = NULL,
## main = NULL, sub = NULL, lty = NULL, lwd = NULL, col = NULL,
## ylim = NULL, ylab = NULL, xlab = NULL, nc, mar.multi = c(0,
## 4, 0, 4), oma.multi = c(6, 4, 6, 4), adj.mtext = NA,
## padj.mtext = NA, col.mtext = NA, ...)
## NULLplot(irf.svara) #fig 2-11
2.2.11 SVAR: Forecast error variance decomposition
## FEVD analysis of SVAR B-type model
args(vars:::fevd.svarest)## function (x, n.ahead = 10, ...)
## NULLfevd.svarb <- fevd(svar.B, n.ahead = 5)
class(fevd.svarb)## [1] "varfevd"methods(class = "varfevd")## [1] plot print
## see '?methods' for accessing help and source codeplot(fevd.svarb)
#add
print(fevd.svarb) #table 2.8## $y1
## y1 y2
## [1,] 1.0000000 0.00000000
## [2,] 0.9803174 0.01968259
## [3,] 0.7538632 0.24613685
## [4,] 0.7538440 0.24615604
## [5,] 0.7459403 0.25405974
##
## $y2
## y1 y2
## [1,] 0.4159728 0.5840272
## [2,] 0.4020796 0.5979204
## [3,] 0.4384861 0.5615139
## [4,] 0.4342112 0.5657888
## [5,] 0.4350269 0.56497312.3 Non-stationary
2.3.1 Stochastic and deterministic trends
set.seed(123456)
e <- rnorm(500)
## pure random walk
rw.nd <- cumsum(e)
## trend
trd <- 1:500
## random walk with drift
rw.wd <- 0.5*trd + cumsum(e)
## deterministic trend and noise
dt <- e + 0.5*trd
## plotting
par(mar=rep(5,4))
plot.ts(dt, lty=1, ylab='', xlab='')
lines(rw.wd, lty=2)
par(new=T)
plot.ts(rw.nd, lty=3, axes=FALSE) #fig 3.1
axis(4, pretty(range(rw.nd)))
lines(rw.nd, lty=3)
legend(10, 18.7, legend=c('det. trend + noise (ls)',
'rw drift (ls)', 'rw (rs)'),
lty=c(1, 2, 3)) 
2.3.2 ARMA versus ARFIMA model
library(fracdiff)## Warning: package 'fracdiff' was built under R version 3.6.3set.seed(123456)
# ARFIMA(0.4,0.4,0.0)
y1 <- fracdiff.sim(n=1000, ar=0.4, ma=0.0, d=0.4)
# ARIMA(0.4,0.0,0.0)
y2 <- arima.sim(model=list(ar=0.4), n=1000)
# Graphics
op <- par(no.readonly=TRUE)
layout(matrix(1:6, 3, 2, byrow=FALSE))
plot.ts(y1$series,
main='Time series plot of long memory',
ylab='')
acf(y1$series, lag.max=100, #fig 3.4
main='Autocorrelations of long memory')
spectrum(y1$series, #fig 3.4
main='Spectral density of long memory')
plot.ts(y2,
main='Time series plot of short memory', ylab='')
acf(y2, lag.max=100, #fig 3.4
main='Autocorrelations of short memory')
spectrum(y2, main='Spectral density of short memory') #fig 3.4
par(op)2.3.3 R/S statistic
library(fracdiff)
set.seed(123456)
# ARFIMA(0.0,0.3,0.0)
y <- fracdiff.sim(n=1000, ar=0.0, ma=0.0, d=0.3)## Warning in min(Mod(polyroot(c(1, -ar)))): no non-missing arguments to min;
## returning Inf# Get the data series, demean this if necessary
y.dm <- y$series
max.y <- max(cumsum(y.dm))
min.y <- min(cumsum(y.dm))
sd.y <- sd(y$series)
RS <- (max.y - min.y)/sd.y
H <- log(RS)/log(1000)
d <- H - 0.52.3.4 Geweke and Porter-Hudak method
library(fracdiff)
set.seed(123456)
y <- fracdiff.sim(n=1000, ar=0.0, ma=0.0, d=0.3)## Warning in min(Mod(polyroot(c(1, -ar)))): no non-missing arguments to min;
## returning Infy.spec <- spectrum(y$series, plot=FALSE)
lhs <- log(y.spec$spec)
rhs <- log(4*(sin(y.spec$freq/2))^2)
gph.reg <- lm(lhs ~ rhs)
gph.sum <- summary(gph.reg)
sqrt(gph.sum$cov.unscaled*pi/6)[2,2]## [1] 0.01663586#add
print(sqrt(gph.sum$cov.unscaled*pi/6)[2,2]) #table 3.1## [1] 0.016635862.4 Cointegration
2.4.1 Spurious regression
library(lmtest)
set.seed(123456)
e1 <- rnorm(500)
e2 <- rnorm(500)
trd <- 1:500
y1 <- 0.8*trd + cumsum(e1)
y2 <- 0.6*trd + cumsum(e2)
sr.reg <- lm(y1 ~ y2)
sr.dw <- dwtest(sr.reg)$statistic
#add
print(summary(sr.reg)) #table 4.1##
## Call:
## lm(formula = y1 ~ y2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.6541 -11.5262 0.3589 11.1423 31.0058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.326969 1.367155 -21.45 <2e-16 ***
## y2 1.440787 0.007519 191.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.71 on 498 degrees of freedom
## Multiple R-squared: 0.9866, Adjusted R-squared: 0.9866
## F-statistic: 3.672e+04 on 1 and 498 DF, p-value: < 2.2e-162.4.2 Engle-Granger procedure with generated data
set.seed(123456)
e1 <- rnorm(100)
e2 <- rnorm(100)
y1 <- cumsum(e1)
y2 <- 0.6*y1 + e2
lr.reg <- lm(y2 ~ y1)
error <- residuals(lr.reg)
error.lagged <- error[-c(1, 100)]
dy1 <- diff(y1)
dy2 <- diff(y2)
diff.dat <- data.frame(embed(cbind(dy1, dy2), 2))
colnames(diff.dat) <- c('dy1', 'dy2', 'dy1.1', 'dy2.1')
ecm.reg <- lm(dy2 ~ error.lagged + dy1.1 + dy2.1,
data=diff.dat)
#add
print(summary(ecm.reg)) #table 4.2 rows 1,2 identical but rows 4,5 are very different - DK why##
## Call:
## lm(formula = dy2 ~ error.lagged + dy1.1 + dy2.1, data = diff.dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9588 -0.5439 0.1370 0.7114 2.3065
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.003398 0.103611 0.033 0.9739
## error.lagged -0.968796 0.158554 -6.110 2.24e-08 ***
## dy1.1 0.245649 0.126996 1.934 0.0561 .
## dy2.1 -0.090117 0.105938 -0.851 0.3971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.026 on 94 degrees of freedom
## Multiple R-squared: 0.5464, Adjusted R-squared: 0.5319
## F-statistic: 37.74 on 3 and 94 DF, p-value: 4.243e-162.4.3 Johansen method with artificially generated data
library(urca)
set.seed(12345)
e1 <- rnorm(250, 0, 0.5)
e2 <- rnorm(250, 0, 0.5)
e3 <- rnorm(250, 0, 0.5)
u1.ar1 <- arima.sim(model = list(ar = 0.75),
innov = e1, n = 250)
u2.ar1 <- arima.sim(model = list(ar = 0.3),
innov = e2, n = 250)
y3 <- cumsum(e3)
y1 <- 0.8 * y3 + u1.ar1
y2 <- -0.3 * y3 + u2.ar1
y.mat <- data.frame(y1, y2, y3)
vecm <- ca.jo(y.mat)
jo.results <- summary(vecm)
vecm.r2 <- cajorls(vecm, r = 2)
class(jo.results)## [1] "sumurca"
## attr(,"package")
## [1] "urca"slotNames(jo.results)## [1] "classname" "test.name" "testreg" "teststat" "cval"
## [6] "bpoint" "signif" "model" "type" "auxstat"
## [11] "lag" "H" "A" "lambda" "pval"
## [16] "V" "W" "P"#add
print(jo.results) #table 4.3, table 4.4##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.27036254 0.15474942 0.01884032
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 2 | 4.72 6.50 8.18 11.65
## r <= 1 | 41.69 12.91 14.90 19.19
## r = 0 | 78.17 18.90 21.07 25.75
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## y1.l2 y2.l2 y3.l2
## y1.l2 1.000000 1.0000000 1.0000000
## y2.l2 -4.732436 0.2273774 0.1513858
## y3.l2 -2.129850 -0.6657324 2.3153224
##
## Weights W:
## (This is the loading matrix)
##
## y1.l2 y2.l2 y3.l2
## y1.d -0.034235501 -0.29705774 -0.008294582
## y2.d 0.145988517 -0.08137695 0.003335110
## y3.d 0.002429191 0.01025326 -0.010523045print(vecm.r2) #table 4.5## $rlm
##
## Call:
## lm(formula = substitute(form1), data = data.mat)
##
## Coefficients:
## y1.d y2.d y3.d
## ect1 -0.331293 0.064612 0.012682
## ect2 0.094473 -0.709385 -0.009165
## constant 0.168371 -0.027019 0.025255
## y1.dl1 -0.227677 0.027012 0.068158
## y2.dl1 0.144452 -0.715607 0.040487
## y3.dl1 0.123467 -0.290828 -0.075251
##
##
## $beta
## ect1 ect2
## y1.l2 1.0000000 0.0000000
## y2.l2 0.0000000 1.0000000
## y3.l2 -0.7328534 0.29519622.4.4 VECM as VAR in levels
library(vars)
vecm.level <- vec2var(vecm, r = 2)
arch.test(vecm.level)##
## ARCH (multivariate)
##
## data: Residuals of VAR object vecm.level
## Chi-squared = 173.46, df = 180, p-value = 0.6231normality.test(vecm.level)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object vecm.level
## Chi-squared = 9.1317, df = 6, p-value = 0.1663
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object vecm.level
## Chi-squared = 5.8516, df = 3, p-value = 0.1191
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object vecm.level
## Chi-squared = 3.2801, df = 3, p-value = 0.3504serial.test(vecm.level)##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object vecm.level
## Chi-squared = 120.92, df = 129, p-value = 0.6817predict(vecm.level)## $y1
## fcst lower upper CI
## [1,] 5.093223 3.817836 6.368610 1.275387
## [2,] 4.867437 3.214760 6.520114 1.652677
## [3,] 4.729502 2.866439 6.592565 1.863063
## [4,] 4.650163 2.631769 6.668558 2.018395
## [5,] 4.607064 2.456642 6.757486 2.150422
## [6,] 4.586909 2.316421 6.857398 2.270489
## [7,] 4.581503 2.198386 6.964620 2.383117
## [8,] 4.585629 2.095338 7.075921 2.490292
## [9,] 4.595931 2.002916 7.188946 2.593015
## [10,] 4.610236 1.918350 7.302123 2.691887
##
## $y2
## fcst lower upper CI
## [1,] -1.575296 -2.528503 -0.6220879 0.9532077
## [2,] -1.542412 -2.596605 -0.4882182 1.0541934
## [3,] -1.555810 -2.658637 -0.4529837 1.1028265
## [4,] -1.577184 -2.721020 -0.4333491 1.1438353
## [5,] -1.596862 -2.778728 -0.4149947 1.1818669
## [6,] -1.613372 -2.831397 -0.3953462 1.2180254
## [7,] -1.627313 -2.880134 -0.3744922 1.2528208
## [8,] -1.639466 -2.925984 -0.3529479 1.2865179
## [9,] -1.650429 -2.969701 -0.3311569 1.3192719
## [10,] -1.660613 -3.011799 -0.3094263 1.3511865
##
## $y3
## fcst lower upper CI
## [1,] 5.430369 4.405761 6.454977 1.024608
## [2,] 5.461409 4.036356 6.886462 1.425053
## [3,] 5.484139 3.752362 7.215916 1.731777
## [4,] 5.508071 3.515104 7.501038 1.992967
## [5,] 5.533684 3.309691 7.757678 2.223993
## [6,] 5.560612 3.127315 7.993909 2.433297
## [7,] 5.588460 2.962420 8.214499 2.626039
## [8,] 5.616920 2.811307 8.422533 2.805613
## [9,] 5.645781 2.671390 8.620171 2.974390
## [10,] 5.674903 2.540795 8.809011 3.134108irf(vecm.level, boot = FALSE)##
## Impulse response coefficients
## $y1
## y1 y2 y3
## [1,] 0.6507197 -0.07111340 0.3165705
## [2,] 0.5313792 -0.09471443 0.3342210
## [3,] 0.4207115 -0.07497365 0.3289089
## [4,] 0.3535795 -0.07483662 0.3261616
## [5,] 0.3111109 -0.08003123 0.3238743
## [6,] 0.2838233 -0.08460734 0.3221977
## [7,] 0.2662038 -0.08786532 0.3210714
## [8,] 0.2547955 -0.09005184 0.3203321
## [9,] 0.2474001 -0.09148793 0.3198502
## [10,] 0.2426041 -0.09242371 0.3195371
## [11,] 0.2394933 -0.09303174 0.3193338
##
## $y2
## y1 y2 y3
## [1,] 0.00000000 0.48111213 -0.1311928
## [2,] 0.05329952 0.17497937 -0.1018415
## [3,] 0.01050817 0.07934508 -0.1156468
## [4,] -0.02396138 0.05610553 -0.1211021
## [5,] -0.04743229 0.04747744 -0.1231880
## [6,] -0.06313329 0.04336982 -0.1243450
## [7,] -0.07343962 0.04110675 -0.1250529
## [8,] -0.08015189 0.03973529 -0.1254997
## [9,] -0.08451240 0.03886819 -0.1257866
## [10,] -0.08734254 0.03831112 -0.1259721
## [11,] -0.08917877 0.03795105 -0.1260922
##
## $y3
## y1 y2 y3
## [1,] 0.00000000 0.0000000 0.3947897
## [2,] 0.04874346 -0.1148161 0.3650815
## [3,] 0.12425300 -0.1240618 0.3612534
## [4,] 0.17707693 -0.1234321 0.3636032
## [5,] 0.21032969 -0.1198584 0.3654301
## [6,] 0.23159439 -0.1164050 0.3667176
## [7,] 0.24532023 -0.1138881 0.3675908
## [8,] 0.25420595 -0.1121907 0.3681660
## [9,] 0.25996535 -0.1110737 0.3685411
## [10,] 0.26370026 -0.1103453 0.3687849
## [11,] 0.26612277 -0.1098718 0.3689432fevd(vecm.level)## $y1
## y1 y2 y3
## [1,] 1.0000000 0.000000000 0.000000000
## [2,] 0.9926630 0.003995459 0.003341587
## [3,] 0.9770177 0.003266243 0.019716067
## [4,] 0.9503105 0.003324248 0.046365265
## [5,] 0.9176063 0.004797534 0.077596120
## [6,] 0.8831519 0.007273673 0.109574452
## [7,] 0.8495802 0.010250495 0.140169301
## [8,] 0.8182411 0.013366630 0.168392272
## [9,] 0.7896646 0.016409187 0.193926229
## [10,] 0.7639236 0.019270134 0.216806301
##
## $y2
## y1 y2 y3
## [1,] 0.02138080 0.9786192 0.00000000
## [2,] 0.04848970 0.9059422 0.04556809
## [3,] 0.06206145 0.8476872 0.09025132
## [4,] 0.07413479 0.7972365 0.12862866
## [5,] 0.08705517 0.7529523 0.15999248
## [6,] 0.10049852 0.7137817 0.18571977
## [7,] 0.11388889 0.6788193 0.20729180
## [8,] 0.12682226 0.6473896 0.22578815
## [9,] 0.13907687 0.6189771 0.24194602
## [10,] 0.15055806 0.5931707 0.25627128
##
## $y3
## y1 y2 y3
## [1,] 0.3667088 0.06297966 0.5703115
## [2,] 0.4008733 0.05217700 0.5469497
## [3,] 0.4100153 0.05246199 0.5375227
## [4,] 0.4124750 0.05379617 0.5337288
## [5,] 0.4126980 0.05498608 0.5323159
## [6,] 0.4121058 0.05596491 0.5319293
## [7,] 0.4112561 0.05676240 0.5319815
## [8,] 0.4103733 0.05741523 0.5322115
## [9,] 0.4095440 0.05795440 0.5325016
## [10,] 0.4087970 0.05840414 0.5327988class(vecm.level)## [1] "vec2var"methods(class = "vec2var")## [1] fevd fitted irf logLik Phi plot predict
## [8] print Psi residuals
## see '?methods' for accessing help and source code#add
print(vecm.level) #should be table 4.6 but deterministic matches, 2 other coeffs exact, rest is different - because reported as A1 A2?##
## Coefficient matrix of lagged endogenous variables:
##
## A1:
## y1.l1 y2.l1 y3.l1
## y1 0.77232286 0.14445173 0.1234669
## y2 0.02701223 0.28439273 -0.2908284
## y3 0.06815811 0.04048737 0.9247494
##
##
## A2:
## y1.l2 y2.l2 y3.l2
## y1 -0.10361609 -0.049978651 0.14721057
## y2 0.03759934 0.006222713 0.03407002
## y3 -0.05547566 -0.049652006 0.06325084
##
##
## Coefficient matrix of deterministic regressor(s).
##
## constant
## y1 0.16837095
## y2 -0.02701889
## y3 0.025255182.5 Order of Integration
2.5.1 ADF test: Integration order for consumption in the United Kingdom
library(urca)
data(Raotbl3)
attach(Raotbl3)
lc <- ts(lc, start=c(1966,4), end=c(1991,2), frequency=4)
lc.ct <- ur.df(lc, lags=3, type='trend')
plot(lc.ct) #fig 5.1
lc.co <- ur.df(lc, lags=3, type='drift')
lc2 <- diff(lc)
lc2.ct <- ur.df(lc2, type='trend', lags=3)
#add
print(summary(lc.ct)) #tables 5.1, 5.2##
## ###############################################
## # 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.044714 -0.006525 0.000129 0.006225 0.045353
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7976591 0.3547775 2.248 0.0270 *
## z.lag.1 -0.0758706 0.0338880 -2.239 0.0277 *
## tt 0.0004915 0.0002159 2.277 0.0252 *
## z.diff.lag1 -0.1063957 0.1006744 -1.057 0.2934
## z.diff.lag2 0.2011373 0.1012373 1.987 0.0500 .
## z.diff.lag3 0.2998586 0.1020548 2.938 0.0042 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01307 on 89 degrees of freedom
## Multiple R-squared: 0.1472, Adjusted R-squared: 0.09924
## F-statistic: 3.071 on 5 and 89 DF, p-value: 0.01325
##
##
## Value of test-statistic is: -2.2389 3.7382 2.5972
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.47print(summary(lc.co)) #tables 5.3, 5.4##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047547 -0.007071 0.000265 0.007731 0.046880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0123237 0.0851358 0.145 0.8852
## z.lag.1 -0.0007356 0.0079043 -0.093 0.9261
## z.diff.lag1 -0.1433015 0.1016454 -1.410 0.1620
## z.diff.lag2 0.1615256 0.1020242 1.583 0.1169
## z.diff.lag3 0.2585280 0.1027364 2.516 0.0136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01337 on 90 degrees of freedom
## Multiple R-squared: 0.09747, Adjusted R-squared: 0.05735
## F-statistic: 2.43 on 4 and 90 DF, p-value: 0.05335
##
##
## Value of test-statistic is: -0.0931 2.8806
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86print(summary(lc2.ct)) #table 5.5 I(2)##
## ###############################################
## # 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.045039 -0.007870 0.000013 0.007807 0.046403
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.864e-03 3.051e-03 1.266 0.2087
## z.lag.1 -8.826e-01 2.013e-01 -4.385 3.2e-05 ***
## tt 3.186e-05 5.112e-05 0.623 0.5348
## z.diff.lag1 -2.253e-01 1.873e-01 -1.203 0.2321
## z.diff.lag2 -4.668e-02 1.600e-01 -0.292 0.7711
## z.diff.lag3 1.775e-01 1.057e-01 1.679 0.0967 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01329 on 88 degrees of freedom
## Multiple R-squared: 0.6147, Adjusted R-squared: 0.5929
## F-statistic: 28.08 on 5 and 88 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -4.3853 6.4477 9.6164
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.472.5.2 PP test: Integration order for consumption in the United Kingdom
library(urca)
data(Raotbl3)
attach(Raotbl3)## The following object is masked _by_ .GlobalEnv:
##
## lc## The following objects are masked from Raotbl3 (pos = 3):
##
## dd682, dd792, dd883, lc, li, lwlc <- ts(lc, start=c(1966,4), end=c(1991,2),
frequency=4)
lc.ct <- ur.pp(lc, type='Z-tau', model='trend',
lags='long')
lc.co <- ur.pp(lc, type='Z-tau', model='constant',
lags='long')
lc2 <- diff(lc)
lc2.ct <- ur.pp(lc2, type='Z-tau', model='trend',
lags='long')
print(summary(lc.ct)) #tables 5.6, 5.7 'constant and trend'##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept and trend
##
##
## Call:
## lm(formula = y ~ y.l1 + trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.045860 -0.006687 0.001637 0.007379 0.047093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5792004 0.3621778 1.599 0.113
## y.l1 0.9468947 0.0335844 28.194 <2e-16 ***
## trend 0.0003472 0.0002155 1.611 0.111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01362 on 95 degrees of freedom
## Multiple R-squared: 0.9946, Adjusted R-squared: 0.9945
## F-statistic: 8695 on 2 and 95 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic, type: Z-tau is: -1.9212
##
## aux. Z statistics
## Z-tau-mu 0.7152
## Z-tau-beta 2.5699
##
## Critical values for Z statistics:
## 1pct 5pct 10pct
## critical values -4.053974 -3.455671 -3.153363print(summary(lc.co)) #tables 5.8, 5.9 'constant only'##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept
##
##
## Call:
## lm(formula = y ~ y.l1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047503 -0.007042 0.001483 0.007879 0.047975
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.010867 0.082462 0.132 0.895
## y.l1 0.999597 0.007643 130.779 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01374 on 96 degrees of freedom
## Multiple R-squared: 0.9944, Adjusted R-squared: 0.9944
## F-statistic: 1.71e+04 on 1 and 96 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic, type: Z-tau is: -0.1285
##
## aux. Z statistics
## Z-tau-mu 0.1993
##
## Critical values for Z statistics:
## 1pct 5pct 10pct
## critical values -3.49776 -2.890909 -2.58224print(summary(lc2.ct)) #table 5.10 'I2'##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept and trend
##
##
## Call:
## lm(formula = y ~ y.l1 + trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.042937 -0.007535 0.000522 0.008133 0.048453
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.340e-03 1.553e-03 4.728 7.96e-06 ***
## y.l1 -1.253e-01 1.025e-01 -1.222 0.225
## trend 1.802e-05 4.996e-05 0.361 0.719
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01377 on 94 degrees of freedom
## Multiple R-squared: 0.01661, Adjusted R-squared: -0.004312
## F-statistic: 0.7939 on 2 and 94 DF, p-value: 0.4551
##
##
## Value of test-statistic, type: Z-tau is: -10.9606
##
## aux. Z statistics
## Z-tau-mu 2.9528
## Z-tau-beta 0.3752
##
## Critical values for Z statistics:
## 1pct 5pct 10pct
## critical values -4.054955 -3.456134 -3.1536332.5.3 ERS tests: Integration order for real GNP in the United States
library(urca)
data(nporg)
gnp <- log(na.omit(nporg[, "gnp.r"]))
gnp.d <- diff(gnp)
gnp.ct.df <- ur.ers(gnp, type = "DF-GLS",
model = "trend", lag.max = 4)
gnp.ct.pt <- ur.ers(gnp, type = "P-test",
model = "trend")
gnp.d.ct.df <- ur.ers(gnp.d, type = "DF-GLS",
model = "trend", lag.max = 4)
gnp.d.ct.pt <- ur.ers(gnp.d, type = "P-test",
model = "trend")
#add
print(summary(gnp.ct.df)) #table 5.11 row 2##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type DF-GLS
## detrending of series with intercept and trend
##
##
## Call:
## lm(formula = dfgls.form, data = data.dfgls)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.166241 -0.035401 0.000747 0.034208 0.140042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## yd.lag -0.14734 0.07086 -2.079 0.04254 *
## yd.diff.lag1 0.39706 0.13566 2.927 0.00507 **
## yd.diff.lag2 0.06112 0.14288 0.428 0.67056
## yd.diff.lag3 -0.07153 0.14241 -0.502 0.61757
## yd.diff.lag4 -0.03216 0.13910 -0.231 0.81806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06071 on 52 degrees of freedom
## Multiple R-squared: 0.2319, Adjusted R-squared: 0.158
## F-statistic: 3.14 on 5 and 52 DF, p-value: 0.01498
##
##
## Value of test-statistic is: -2.0793
##
## Critical values of DF-GLS are:
## 1pct 5pct 10pct
## critical values -3.58 -3.03 -2.74print(summary(gnp.ct.pt)) #table 5.11 row 1##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type P-test
## detrending of series with intercept and trend
##
## Value of test-statistic is: 6.6528
##
## Critical values of P-test are:
## 1pct 5pct 10pct
## critical values 4.26 5.64 6.79print(summary(gnp.d.ct.df)) #table 5.12 row 2##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type DF-GLS
## detrending of series with intercept and trend
##
##
## Call:
## lm(formula = dfgls.form, data = data.dfgls)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.155873 -0.026967 0.000705 0.043680 0.152529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## yd.lag -1.04464 0.25305 -4.128 0.000135 ***
## yd.diff.lag1 0.35713 0.21614 1.652 0.104619
## yd.diff.lag2 0.33770 0.18699 1.806 0.076826 .
## yd.diff.lag3 0.18645 0.16535 1.128 0.264758
## yd.diff.lag4 0.09155 0.13831 0.662 0.511017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06277 on 51 degrees of freedom
## Multiple R-squared: 0.3795, Adjusted R-squared: 0.3187
## F-statistic: 6.239 on 5 and 51 DF, p-value: 0.0001381
##
##
## Value of test-statistic is: -4.1283
##
## Critical values of DF-GLS are:
## 1pct 5pct 10pct
## critical values -3.58 -3.03 -2.74print(summary(gnp.d.ct.pt)) #table 5.12 row 1##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type P-test
## detrending of series with intercept and trend
##
## Value of test-statistic is: 2.6761
##
## Critical values of P-test are:
## 1pct 5pct 10pct
## critical values 4.26 5.64 6.792.5.4 SP test: Integration order for nominal GNP of the United States
library(urca)
data(nporg)
gnp <- na.omit(nporg[, "gnp.n"])
gnp.tau.sp <- ur.sp(gnp, type = "tau", pol.deg=2,
signif=0.05)
gnp.rho.sp <- ur.sp(gnp, type = "rho", pol.deg=2,
signif=0.05)
#--added
plot(gnp,ty='l') #fig 5.2
print(summary(gnp.tau.sp)) #tables 5.13, 5.14 row 1##
## ###################################
## # Schmidt-Phillips Unit Root Test #
## ###################################
##
##
## Call:
## lm(formula = sp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24494 -5959 1334 7556 23919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9355.40157 7395.30256 1.265 0.2110
## y.lagged 0.98813 0.04112 24.033 <2e-16 ***
## trend.exp1 -982.32296 610.13476 -1.610 0.1129
## trend.exp2 30.27619 16.15727 1.874 0.0661 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12220 on 57 degrees of freedom
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.9977
## F-statistic: 8538 on 3 and 57 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -2.0257
## Critical value for a significance level of 0.05
## is: -3.65print(summary(gnp.rho.sp)) #tables 5.13, 5.14 row 2##
## ###################################
## # Schmidt-Phillips Unit Root Test #
## ###################################
##
##
## Call:
## lm(formula = sp.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24494 -5959 1334 7556 23919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9355.40157 7395.30256 1.265 0.2110
## y.lagged 0.98813 0.04112 24.033 <2e-16 ***
## trend.exp1 -982.32296 610.13476 -1.610 0.1129
## trend.exp2 30.27619 16.15727 1.874 0.0661 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12220 on 57 degrees of freedom
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.9977
## F-statistic: 8538 on 3 and 57 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -8.5056
## Critical value for a significance level of 0.05
## is: -23.72.5.5 KPSS test: Integration order for interest rates and nominal wages in the United States
library(urca)
data(nporg)
ir <- na.omit(nporg[, "bnd"])
wg <- log(na.omit(nporg[, "wg.n"]))
ir.kpss <- ur.kpss(ir, type = "mu", use.lag=8)
wg.kpss <- ur.kpss(wg, type = "tau", use.lag=8)
#add
print(summary(ir.kpss)) #table 5.15 row 1##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 8 lags.
##
## Value of test-statistic is: 0.1325
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739print(summary(wg.kpss)) #table 5.15 row 2##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 8 lags.
##
## Value of test-statistic is: 0.1007
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.2162.6 Further Tests
2.6.1 Random walk with drift and structural break
set.seed(123456)
e <- rnorm(500)
## trend
trd <- 1:500
S <- c(rep(0, 249), rep(1, 251))
## random walk with drift
y1 <- 0.1*trd + cumsum(e)
## random walk with drift and shift
y2 <- 0.1*trd + 10*S + cumsum(e)
#-add
plot(as.zoo(cbind(y1,y2)),scr=1,lty=c(1,2)) #fig 6.1
2.6.2 Unit roots and structural break: Zivot-Andrews test
library(urca)
data(nporg)
wg.n <- log(na.omit(nporg[, "wg.n"]))
za.wg.n <- ur.za(wg.n, model = "intercept", lag = 7)
plot(za.wg.n)
wg.r <- log(na.omit(nporg[, "wg.r"])) #fig 6.2
za.wg.r <- ur.za(wg.r, model = "both", lag = 8)
plot(za.wg.r)
#add
print(summary(za.wg.n)) #table 6.1, 6.3 row 1##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.106309 -0.021981 -0.001136 0.031850 0.137894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.987820 0.372428 5.337 2.09e-06 ***
## y.l1 0.660045 0.064113 10.295 3.77e-14 ***
## trend 0.017321 0.003256 5.319 2.23e-06 ***
## y.dl1 0.497942 0.112122 4.441 4.70e-05 ***
## y.dl2 0.055750 0.130817 0.426 0.67174
## y.dl3 0.149426 0.127817 1.169 0.24771
## y.dl4 0.061080 0.126575 0.483 0.63144
## y.dl5 0.006119 0.126417 0.048 0.96158
## y.dl6 0.141893 0.124862 1.136 0.26100
## y.dl7 0.267111 0.119471 2.236 0.02968 *
## du -0.160825 0.038681 -4.158 0.00012 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05367 on 52 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.996, Adjusted R-squared: 0.9953
## F-statistic: 1305 on 10 and 52 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -5.3024
## Critical values: 0.01= -5.34 0.05= -4.8 0.1= -4.58
##
## Potential break point at position: 30print(summary(za.wg.r)) #table 6.2, 6.3 row 2##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.072564 -0.010853 0.002307 0.015563 0.058436
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.567141 0.532696 4.819 1.43e-05 ***
## y.l1 0.114607 0.186630 0.614 0.542001
## trend 0.012364 0.002755 4.488 4.37e-05 ***
## y.dl1 0.611066 0.166238 3.676 0.000588 ***
## y.dl2 0.351601 0.168613 2.085 0.042278 *
## y.dl3 0.441312 0.156769 2.815 0.007005 **
## y.dl4 0.256360 0.145265 1.765 0.083835 .
## y.dl5 0.138055 0.134589 1.026 0.310045
## y.dl6 0.059095 0.126183 0.468 0.641627
## y.dl7 0.167333 0.120061 1.394 0.169688
## y.dl8 0.148557 0.121008 1.228 0.225439
## du 0.084907 0.019616 4.328 7.39e-05 ***
## dt 0.008103 0.002201 3.682 0.000577 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02996 on 49 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.9953, Adjusted R-squared: 0.9941
## F-statistic: 860.7 on 12 and 49 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -4.7441
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 412.6.3 HEGY test for seasonal unit roots
#as well as having renamed hegy lowercase, the args have changed a lot so this is not run [not relevant, seasonal]
library(urca)
library(uroot)
data(UKconinc)
incl <- ts(UKconinc$incl, start = c(1955,1),
end = c(1984,4), frequency = 4)
hegy000 <- hegy.test(wts = incl, itsd = c(0, 0, c(0)),
selectlags = list(mode = c(1,4,5)))
hegy100 <- hegy.test(wts = incl, itsd = c(1, 0, c(0)),
selectlags = list(mode = c(1,4,5)))
hegy110 <- hegy.test(wts = incl, itsd = c(1, 1, c(0)),
selectlags = list(mode = c(1,4,5)))
hegy101 <- hegy.test(wts = incl,
itsd = c(1, 0, c(1, 2, 3)),
selectlags = list(mode = c(1,4,5)))
hegy111 <- hegy.test(wts = incl,
itsd = c(1, 1, c(1, 2, 3)),
selectlags = list(mode = c(1,4,5)))2.7 Single-Equation
2.7.1 Engle-Granger: Long-run relationship of consumption, income, and wealth in the United Kingdom
library(tseries)
library(urca)
data(Raotbl3)
attach(Raotbl3)## The following object is masked _by_ .GlobalEnv:
##
## lc## The following objects are masked from Raotbl3 (pos = 4):
##
## dd682, dd792, dd883, lc, li, lw## The following objects are masked from Raotbl3 (pos = 5):
##
## dd682, dd792, dd883, lc, li, lwlc <- ts(lc, start=c(1966,4), end=c(1991,2),
frequency=4)
li <- ts(li, start=c(1966,4), end=c(1991,2),
frequency=4)
lw <- ts(lw, start=c(1966,4), end=c(1991,2),
frequency=4)
ukcons <- window(cbind(lc, li, lw), start=c(1967, 2),
end=c(1991,2))
lc.eq <- summary(lm(lc ~ li + lw, data=ukcons))
li.eq <- summary(lm(li ~ lc + lw, data=ukcons))
lw.eq <- summary(lm(lw ~ li + lc, data=ukcons))
error.lc <- ts(resid(lc.eq), start=c(1967,2),
end=c(1991,2), frequency=4)
error.li <- ts(resid(li.eq), start=c(1967,2),
end=c(1991,2), frequency=4)
error.lw <- ts(resid(lw.eq), start=c(1967,2),
end=c(1991,2), frequency=4)
ci.lc <- ur.df(error.lc, lags=1, type='none')
ci.li <- ur.df(error.li, lags=1, type='none')
ci.lw <- ur.df(error.lw, lags=1, type='none')
jb.lc <- jarque.bera.test(error.lc)
jb.li <- jarque.bera.test(error.li)
jb.lw <- jarque.bera.test(error.lw)
#add
print(summary(ci.lc)@testreg$coefficients[1,3]) #table 7.1 row 1 col 1## [1] -4.141447print(summary(ci.li)@testreg$coefficients[1,3]) #table 7.1 row 2 col 1## [1] -4.059031print(summary(ci.lw)@testreg$coefficients[1,3]) #table 7.1 row 3 col 1## [1] -2.713707print(jb.lc) #table 7.1 row 1 col 2,3##
## Jarque Bera Test
##
## data: error.lc
## X-squared = 0.66447, df = 2, p-value = 0.7173print(jb.li) #table 7.1 row 2 col 2,3##
## Jarque Bera Test
##
## data: error.li
## X-squared = 0.065267, df = 2, p-value = 0.9679print(jb.lw) #table 7.1 row 3 col 2,3##
## Jarque Bera Test
##
## data: error.lw
## X-squared = 3.253, df = 2, p-value = 0.19662.7.2 Engle-Granger: ECM for consumption and income of the United Kingdom
ukcons2 <- ts(embed(diff(ukcons), dim=2),
start=c(1967,4), freq=4)
colnames(ukcons2) <- c('lc.d', 'li.d', 'lw.d',
'lc.d1', 'li.d1', 'lw.d1')
error.ecm1 <- window(lag(error.lc, k=-1),
start=c(1967,4), end=c(1991, 2))
error.ecm2 <- window(lag(error.li, k=-1),
start=c(1967,4), end=c(1991, 2))
ecm.eq1 <- lm(lc.d ~ error.ecm1 + lc.d1 + li.d1 + lw.d1,
data=ukcons2)
ecm.eq2 <- lm(li.d ~ error.ecm2 + lc.d1 + li.d1 + lw.d1,
data=ukcons2)
#add
print(summary(ecm.eq1)) #table 7.2##
## Call:
## lm(formula = lc.d ~ error.ecm1 + lc.d1 + li.d1 + lw.d1, data = ukcons2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.044333 -0.006459 0.000915 0.006185 0.043750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005770 0.001496 3.856 0.000216 ***
## error.ecm1 0.062528 0.098408 0.635 0.526785
## lc.d1 -0.285563 0.115824 -2.465 0.015579 *
## li.d1 0.261417 0.086361 3.027 0.003221 **
## lw.d1 0.082694 0.031688 2.610 0.010614 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01282 on 90 degrees of freedom
## Multiple R-squared: 0.1704, Adjusted R-squared: 0.1335
## F-statistic: 4.621 on 4 and 90 DF, p-value: 0.001938print(summary(ecm.eq2)) #table 7.3##
## Call:
## lm(formula = li.d ~ error.ecm2 + lc.d1 + li.d1 + lw.d1, data = ukcons2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.048578 -0.007692 0.000699 0.007625 0.049778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.006643 0.001880 3.535 0.000647 ***
## error.ecm2 -0.539482 0.114211 -4.724 8.49e-06 ***
## lc.d1 -0.149635 0.146448 -1.022 0.309633
## li.d1 -0.006035 0.108544 -0.056 0.955784
## lw.d1 0.062673 0.039784 1.575 0.118694
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0161 on 90 degrees of freedom
## Multiple R-squared: 0.2499, Adjusted R-squared: 0.2166
## F-statistic: 7.497 on 4 and 90 DF, p-value: 2.937e-052.7.3 Phillips-Ouliaris: Long-run relationship of consumption, income, and wealth in the United Kingdom
library(urca)
data(Raotbl3)
attach(Raotbl3)## The following objects are masked _by_ .GlobalEnv:
##
## lc, li, lw## The following objects are masked from Raotbl3 (pos = 3):
##
## dd682, dd792, dd883, lc, li, lw## The following objects are masked from Raotbl3 (pos = 5):
##
## dd682, dd792, dd883, lc, li, lw## The following objects are masked from Raotbl3 (pos = 6):
##
## dd682, dd792, dd883, lc, li, lwlc <- ts(lc, start=c(1966,4), end=c(1991,2),
frequency=4)
li <- ts(li, start=c(1966,4), end=c(1991,2),
frequency=4)
lw <- ts(lw, start=c(1966,4), end=c(1991,2),
frequency=4)
ukcons <- window(cbind(lc, li, lw), start=c(1967, 2),
end=c(1991,2))
pu.test <- summary(ca.po(ukcons, demean='const',
type='Pu'))
pz.test <- summary(ca.po(ukcons, demean='const',
type='Pz'))
#add
print(pu.test) #table 7.4 row 1##
## ########################################
## # Phillips and Ouliaris Unit Root Test #
## ########################################
##
## Test of type Pu
## detrending of series with constant only
##
##
## Call:
## lm(formula = z[, 1] ~ z[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.037809 -0.009914 0.000317 0.007332 0.042950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.178458 0.103491 -1.724 0.0879 .
## z[, -1]li 0.910971 0.013457 67.697 <2e-16 ***
## z[, -1]lw 0.079761 0.007764 10.274 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01656 on 94 degrees of freedom
## Multiple R-squared: 0.9918, Adjusted R-squared: 0.9916
## F-statistic: 5700 on 2 and 94 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: 58.9108
##
## Critical values of Pu are:
## 10pct 5pct 1pct
## critical values 33.6955 40.5252 53.8731print(pz.test) #table 7.4 row 2##
## ########################################
## # Phillips and Ouliaris Unit Root Test #
## ########################################
##
## Test of type Pz
## detrending of series with constant only
##
## Response lc :
##
## Call:
## lm(formula = lc ~ zr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047350 -0.007268 0.001987 0.007124 0.049443
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.022293 0.087754 -0.254 0.80003
## zrlc 0.853380 0.085127 10.025 < 2e-16 ***
## zrli 0.118152 0.078442 1.506 0.13544
## zrlw 0.024679 0.009372 2.633 0.00992 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01351 on 92 degrees of freedom
## Multiple R-squared: 0.9945, Adjusted R-squared: 0.9943
## F-statistic: 5566 on 3 and 92 DF, p-value: < 2.2e-16
##
##
## Response li :
##
## Call:
## lm(formula = li ~ zr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043415 -0.009629 -0.000524 0.008431 0.051594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.13002 0.10555 1.232 0.2211
## zrlc 0.49383 0.10239 4.823 5.58e-06 ***
## zrli 0.52867 0.09435 5.603 2.18e-07 ***
## zrlw -0.02429 0.01127 -2.155 0.0338 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01625 on 92 degrees of freedom
## Multiple R-squared: 0.9918, Adjusted R-squared: 0.9915
## F-statistic: 3697 on 3 and 92 DF, p-value: < 2.2e-16
##
##
## Response lw :
##
## Call:
## lm(formula = lw ~ zr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.162475 -0.016750 0.004733 0.028571 0.086684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.31857 0.27819 -1.145 0.255
## zrlc 0.13079 0.26986 0.485 0.629
## zrli -0.07984 0.24867 -0.321 0.749
## zrlw 0.98357 0.02971 33.106 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04283 on 92 degrees of freedom
## Multiple R-squared: 0.9815, Adjusted R-squared: 0.9809
## F-statistic: 1627 on 3 and 92 DF, p-value: < 2.2e-16
##
##
##
## Value of test-statistic is: 88.0345
##
## Critical values of Pz are:
## 10pct 5pct 1pct
## critical values 80.2034 89.7619 109.45252.8 Multiple-Equation
2.8.1 Johansen-Juselius Unrestricted Cointegration
library(urca)
data(UKpppuip)
names(UKpppuip)## [1] "p1" "p2" "e12" "i1" "i2" "doilp0" "doilp1"attach(UKpppuip)
dat1 <- cbind(p1, p2, e12, i1, i2)
dat2 <- cbind(doilp0, doilp1)
args('ca.jo')## function (x, type = c("eigen", "trace"), ecdet = c("none", "const",
## "trend"), K = 2, spec = c("longrun", "transitory"), season = NULL,
## dumvar = NULL)
## NULLH1 <- ca.jo(dat1, type = 'trace', K = 2, season = 4,
dumvar = dat2)
H1.trace <- summary(ca.jo(dat1, type = 'trace', K = 2,
season = 4, dumvar = dat2))
H1.eigen <- summary(ca.jo(dat1, type = 'eigen', K = 2,
season = 4, dumvar = dat2))
#-add
plot(as.zoo(cbind(dat1,dat2[,1,drop=F]))) #fig 8.1
print(H1.eigen) #table 8.1, 8.3 'eigenvectors'##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.40672818 0.28538240 0.25415335 0.10230406 0.08287097
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 4 | 5.19 6.50 8.18 11.65
## r <= 3 | 6.48 12.91 14.90 19.19
## r <= 2 | 17.59 18.90 21.07 25.75
## r <= 1 | 20.16 24.78 27.14 32.14
## r = 0 | 31.33 30.84 33.32 38.78
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## p1.l2 p2.l2 e12.l2 i1.l2 i2.l2
## p1.l2 1.0000000 1.000000 1.000000 1.0000000 1.0000000
## p2.l2 -0.9086265 -1.143047 -1.272628 -2.4001444 -1.4528820
## e12.l2 -0.9321133 -3.363042 1.113631 1.1221619 -0.4805235
## i1.l2 -3.3746393 35.243576 2.746828 -0.4088865 2.2775510
## i2.l2 -1.8906210 -32.917370 -2.835714 2.9863624 0.7628011
##
## Weights W:
## (This is the loading matrix)
##
## p1.l2 p2.l2 e12.l2 i1.l2 i2.l2
## p1.d -0.06816507 0.0011795779 -0.002790218 0.001373599 -0.01333013
## p2.d -0.01773477 0.0001220008 -0.014159241 0.013178503 0.00755575
## e12.d 0.10065321 -0.0001432122 -0.055628059 -0.035400025 -0.04707585
## i1.d 0.03434737 -0.0041631581 -0.010363374 0.012309982 -0.02394672
## i2.d 0.05766426 0.0082830953 0.004821036 0.026984801 -0.01006765print(H1.trace) #table 8.2, 8.4 'weights'##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.40672818 0.28538240 0.25415335 0.10230406 0.08287097
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 4 | 5.19 6.50 8.18 11.65
## r <= 3 | 11.67 15.66 17.95 23.52
## r <= 2 | 29.26 28.71 31.52 37.22
## r <= 1 | 49.42 45.23 48.28 55.43
## r = 0 | 80.75 66.49 70.60 78.87
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## p1.l2 p2.l2 e12.l2 i1.l2 i2.l2
## p1.l2 1.0000000 1.000000 1.000000 1.0000000 1.0000000
## p2.l2 -0.9086265 -1.143047 -1.272628 -2.4001444 -1.4528820
## e12.l2 -0.9321133 -3.363042 1.113631 1.1221619 -0.4805235
## i1.l2 -3.3746393 35.243576 2.746828 -0.4088865 2.2775510
## i2.l2 -1.8906210 -32.917370 -2.835714 2.9863624 0.7628011
##
## Weights W:
## (This is the loading matrix)
##
## p1.l2 p2.l2 e12.l2 i1.l2 i2.l2
## p1.d -0.06816507 0.0011795779 -0.002790218 0.001373599 -0.01333013
## p2.d -0.01773477 0.0001220008 -0.014159241 0.013178503 0.00755575
## e12.d 0.10065321 -0.0001432122 -0.055628059 -0.035400025 -0.04707585
## i1.d 0.03434737 -0.0041631581 -0.010363374 0.012309982 -0.02394672
## i2.d 0.05766426 0.0082830953 0.004821036 0.026984801 -0.010067652.8.2 H1 model: Transformations and cointegration relations
beta <- H1@V
beta[,2] <- beta[,2]/beta[4,2]
beta[,3] <- beta[,3]/beta[4,3]
alpha <- H1@PI%*%solve(t(beta))
beta1 <- cbind(beta[,1:2], H1@V[,3:5])
ci.1 <- ts((H1@x%*%beta1)[-c(1,2),], start=c(1972, 3),
end=c(1987, 2), frequency=4)
ci.2 <- ts(H1@RK%*%beta1, start=c(1972, 3),
end=c(1987, 2), frequency=4)
#-add
plot(cbind(ci.1[,1:2],ci.2[,1:2]),nc=2) # fig 8.2
print(H1@PI) #resembles table 8.5, labels right, and 'PI' sound correct ... but not the same## p1.l2 p2.l2 e12.l2 i1.l2 i2.l2
## p1.d -0.081732241 0.080209451 0.064410165 0.23301914 0.09189181
## p2.d -0.011037753 -0.008613652 0.011510041 0.03707528 0.11478466
## e12.d -0.037593941 0.132862187 -0.172391224 -0.59025936 -0.16946445
## i1.d 0.008184099 -0.008015502 -0.004234984 -0.35067421 0.11998541
## i2.d 0.087685539 -0.118138928 -0.041118105 0.07660906 -0.322443282.8.3 H4 model: Testing for weak exogenity
A1 <- matrix(c(1,0,0,0,0, 0,0,1,0,0,
0,0,0,1,0, 0,0,0,0,1),
nrow=5, ncol=4)
A2 <- matrix(c(1,0,0,0,0, 0,1,0,0,0,
0,0,1,0,0, 0,0,0,1,0),
nrow=5, ncol=4)
H41 <- summary(alrtest(z = H1, A = A1, r = 2))
H42 <- summary(alrtest(z = H1, A = A2, r = 2))
#add
print(H41) #table 8.6 row 1##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under linear restrictions on beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 0 0 0
## [3,] 0 1 0 0
## [4,] 0 0 1 0
## [5,] 0 0 0 1
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.4002 0.2854 0.1674 0.0880 0.0000
##
## The value of the likelihood ratio test statistic:
## 0.66 distributed as chi square with 2 df.
## The p-value of the test statistic is: 0.72
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2]
## RK.p1.l2 1.0000 1.0000
## RK.p2.l2 -0.9249 -1.1420
## RK.e12.l2 -0.9690 -3.2247
## RK.i1.l2 -3.4812 34.6544
## RK.i2.l2 -1.7651 -32.3316
##
## Weights W of the restricted VAR:
##
## [,1] [,2]
## [1,] -0.0645 0.0012
## [2,] 0.0000 0.0000
## [3,] 0.1259 -0.0005
## [4,] 0.0403 -0.0043
## [5,] 0.0566 0.0085print(H42) #table 8.6 row 2##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under linear restrictions on beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1
## [5,] 0 0 0 0
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.3870 0.2560 0.1945 0.0856 0.0000
##
## The value of the likelihood ratio test statistic:
## 4.38 distributed as chi square with 2 df.
## The p-value of the test statistic is: 0.11
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2]
## RK.p1.l2 1.0000 1.0000
## RK.p2.l2 -0.8518 -1.4443
## RK.e12.l2 -0.9525 2.1238
## RK.i1.l2 -4.6381 -0.7205
## RK.i2.l2 -1.0922 1.3813
##
## Weights W of the restricted VAR:
##
## [,1] [,2]
## [1,] -0.0624 -0.0022
## [2,] -0.0157 -0.0101
## [3,] 0.1029 -0.0418
## [4,] 0.0363 -0.0049
## [5,] 0.0000 0.00002.8.4 H3 model: Testing for restrictions in all cointegration relations
H.31 <- matrix(c(1,-1,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),
c(5,3))
H.32 <- matrix(c(1,0,0,0,0, 0,1,0,0,0, 0,0,1,0,0,
0,0,0,1,-1), c(5,4))
H31 <- summary(blrtest(z = H1, H = H.31, r = 2))
H32 <- summary(blrtest(z = H1, H = H.32, r = 2))
#add
print(H31) #table 8.7 row 1##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under linear restrictions on beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] -1 0 0
## [3,] -1 0 0
## [4,] 0 1 0
## [5,] 0 0 1
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.3855 0.2776 0.0895
##
## The value of the likelihood ratio test statistic:
## 2.76 distributed as chi square with 4 df.
## The p-value of the test statistic is: 0.6
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2] [,3]
## [1,] 1.0000 1.0000 1.0000
## [2,] -1.0000 -1.0000 -1.0000
## [3,] -1.0000 -1.0000 -1.0000
## [4,] -2.6141 28.7802 2.5168
## [5,] -2.0949 -26.2994 -0.1954
##
## Weights W of the restricted VAR:
##
## [,1] [,2] [,3]
## p1.d -0.0577 0.0012 -0.0115
## p2.d 0.0004 -0.0030 -0.0013
## e12.d 0.1323 -0.0121 -0.0353
## i1.d 0.0321 -0.0069 -0.0200
## i2.d 0.0683 0.0094 -0.0136print(H32) #table 8.7 row 2##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under linear restrictions on beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1
## [5,] 0 0 0 -1
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.2857 0.2542 0.1458 0.0927
##
## The value of the likelihood ratio test statistic:
## 13.71 distributed as chi square with 2 df.
## The p-value of the test statistic is: 0
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2] [,3] [,4]
## [1,] 1.0000 1.0000 1.0000 1.0000
## [2,] -1.1268 -1.2754 -1.6080 1.7357
## [3,] -2.3411 1.1611 -0.0469 -5.6116
## [4,] 20.7157 2.9055 -0.7435 7.3346
## [5,] -20.7157 -2.9055 0.7435 -7.3346
##
## Weights W of the restricted VAR:
##
## [,1] [,2] [,3] [,4]
## p1.d 0.0013 -0.0026 -0.0460 -0.0029
## p2.d 0.0001 -0.0138 0.0181 -0.0008
## e12.d 0.0004 -0.0542 -0.0558 -0.0018
## i1.d -0.0065 -0.0101 0.0077 -0.0053
## i2.d 0.0139 0.0046 0.0502 -0.00492.8.5 H3 model: Testing for partly known cointegration relations
H.51 <- c(1, -1, -1, 0, 0)
H.52 <- c(0, 0, 0, 1, -1)
H51 <- summary(bh5lrtest(z = H1, H = H.51, r = 2))
H52 <- summary(bh5lrtest(z = H1, H = H.52, r = 2))
#add
print(H51) #table 8.8 row 1##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under partly known beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1]
## [1,] 1
## [2,] -1
## [3,] -1
## [4,] 0
## [5,] 0
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.3956 0.2812 0.2541 0.1008
##
## The value of the likelihood ratio test statistic:
## 14.52 distributed as chi square with 3 df.
## The p-value of the test statistic is: 0
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2]
## [1,] 1 1.0000
## [2,] -1 0.7845
## [3,] -1 0.2155
## [4,] 0 -64.9725
## [5,] 0 -40.7031
##
## Weights W of the restricted VAR:
##
## [,1] [,2]
## p1.d -0.0744 -0.0031
## p2.d -0.0147 -0.0009
## e12.d 0.0724 0.0057
## i1.d 0.0164 0.0020
## i2.d 0.0416 0.0035print(H52) #table 8.8 row 2##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under partly known beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 1
## [5,] -1
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.4064 0.2607 0.1052 0.1007
##
## The value of the likelihood ratio test statistic:
## 1.89 distributed as chi square with 3 df.
## The p-value of the test statistic is: 0.59
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2]
## [1,] 0 1.0000
## [2,] 0 -0.9102
## [3,] 0 -0.9334
## [4,] 1 -2.6163
## [5,] -1 -2.6163
##
## Weights W of the restricted VAR:
##
## [,1] [,2]
## p1.d 0.0881 -0.0678
## p2.d -0.0597 -0.0189
## e12.d -0.2731 0.0975
## i1.d -0.1450 0.0320
## i2.d 0.2197 0.06252.8.6 H6 model: Testing of restrictions on r1 cointegration relations
H.6 <- matrix(rbind(diag(3), c(0, 0, 0), c(0, 0, 0)),
nrow=5, ncol=3)
H6 <- summary(bh6lrtest(z = H1, H = H.6,
r = 2, r1 = 1))
#add
print(H6) #table 8.9##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under partly known beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
## [4,] 0 0 0
## [5,] 0 0 0
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.4067 0.2807 0.1493 0.0906
##
## The value of the likelihood ratio test statistic:
## 4.93 distributed as chi square with 1 df.
## The p-value of the test statistic is: 0.03
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1] [,2]
## [1,] 1.0000 1.0000
## [2,] -1.2245 -4.4520
## [3,] 0.1303 -3.8087
## [4,] 0.0000 -13.9017
## [5,] 0.0000 -7.7497
##
## Weights W of the restricted VAR:
##
## [,1] [,2]
## p1.d -0.0592 -0.0039
## p2.d -0.0166 0.0006
## e12.d -0.0879 0.0018
## i1.d 0.0041 0.0023
## i2.d 0.0881 0.00862.8.7 H1 model: Inference on cointegration rank for Danish money demand function allowing for structural shift
data(denmark)
sjd <- denmark[, c("LRM", "LRY", "IBO", "IDE")]
sjd.vecm <- summary(ca.jo(sjd, ecdet = "const",
type = "eigen",
K = 2,
spec = "longrun",
season = 4))
lue.vecm <- summary(cajolst(sjd, season=4))
#add
print(sjd.vecm) #table 8.10##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 4.331654e-01 1.775836e-01 1.127905e-01 4.341130e-02 6.927550e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 2.35 7.52 9.24 12.97
## r <= 2 | 6.34 13.75 15.67 20.20
## r <= 1 | 10.36 19.77 22.00 26.81
## r = 0 | 30.09 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## LRM.l2 LRY.l2 IBO.l2 IDE.l2 constant
## LRM.l2 1.000000 1.0000000 1.0000000 1.000000 1.0000000
## LRY.l2 -1.032949 -1.3681031 -3.2266580 -1.883625 -0.6336946
## IBO.l2 5.206919 0.2429825 0.5382847 24.399487 1.6965828
## IDE.l2 -4.215879 6.8411103 -5.6473903 -14.298037 -1.8951589
## constant -6.059932 -4.2708474 7.8963696 -2.263224 -8.0330127
##
## Weights W:
## (This is the loading matrix)
##
## LRM.l2 LRY.l2 IBO.l2 IDE.l2 constant
## LRM.d -0.21295494 -0.00481498 0.035011128 2.028908e-03 -1.726523e-13
## LRY.d 0.11502204 0.01975028 0.049938460 1.108654e-03 9.428195e-14
## IBO.d 0.02317724 -0.01059605 0.003480357 -1.573742e-03 4.143714e-14
## IDE.d 0.02941109 -0.03022917 -0.002811506 -4.767627e-05 7.781415e-14print(lue.vecm) #table 8.11##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in shift correction
##
## Eigenvalues (lambda):
## [1] 0.42098147 0.27098513 0.17330604 0.06127991
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.15 5.42 6.79 10.04
## r <= 2 | 11.62 13.78 15.83 19.85
## r <= 1 | 24.33 25.93 28.45 33.76
## r = 0 | 42.95 42.08 45.20 51.60
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## LRM.l1 LRY.l1 IBO.l1 IDE.l1
## LRM.l1 1.0000000 1.000000 1.0000000 1.0000000
## LRY.l1 0.8490645 -2.319498 -1.8603632 -0.2148260
## IBO.l1 6.0747337 6.824430 0.8571653 0.4685797
## IDE.l1 1.0858006 -11.662840 2.7346259 -0.5253649
##
## Weights W:
## (This is the loading matrix)
##
## LRM.l1 LRY.l1 IBO.l1 IDE.l1
## LRM -0.146180842 -0.034940907 0.04306167 -0.066397980
## LRY -0.081222325 0.090040394 0.10130486 -0.011416802
## IBO -0.021090906 0.002883555 -0.01985965 0.030806000
## IDE 0.004947721 0.022143257 -0.03940262 0.0080278782.8.8 Canadian data set: Preliminary analysis
library(vars)
data(Canada)
summary(Canada)## e prod rw U
## Min. :928.6 Min. :401.3 Min. :386.1 Min. : 6.700
## 1st Qu.:935.4 1st Qu.:404.8 1st Qu.:423.9 1st Qu.: 7.782
## Median :946.0 Median :406.5 Median :444.4 Median : 9.450
## Mean :944.3 Mean :407.8 Mean :440.8 Mean : 9.321
## 3rd Qu.:950.0 3rd Qu.:410.7 3rd Qu.:461.1 3rd Qu.:10.607
## Max. :961.8 Max. :418.0 Max. :470.0 Max. :12.770plot(Canada, nc = 2) #fig 8.3
2.8.9 Canadian data set: ADF-test regressions
summary(ur.df(Canada[, "prod"],
type = "trend", lags = 2))##
## ###############################################
## # 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
## -2.19924 -0.38994 0.04294 0.41914 1.71660
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.415228 15.309403 1.987 0.0506 .
## z.lag.1 -0.075791 0.038134 -1.988 0.0505 .
## tt 0.013896 0.006422 2.164 0.0336 *
## z.diff.lag1 0.284866 0.114359 2.491 0.0149 *
## z.diff.lag2 0.080019 0.116090 0.689 0.4927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6851 on 76 degrees of freedom
## Multiple R-squared: 0.1354, Adjusted R-squared: 0.08993
## F-statistic: 2.976 on 4 and 76 DF, p-value: 0.02438
##
##
## Value of test-statistic is: -1.9875 2.3 2.3817
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.47summary(ur.df(diff(Canada[, "prod"]),
type = "drift", lags = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.05124 -0.39530 0.07819 0.41109 1.75129
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.11534 0.08029 1.437 0.155
## z.lag.1 -0.68893 0.13350 -5.160 1.83e-06 ***
## z.diff.lag -0.04274 0.11275 -0.379 0.706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6971 on 78 degrees of freedom
## Multiple R-squared: 0.3615, Adjusted R-squared: 0.3451
## F-statistic: 22.08 on 2 and 78 DF, p-value: 2.526e-08
##
##
## Value of test-statistic is: -5.1604 13.3184
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86summary(ur.df(Canada[, "e"],
type = "trend", lags = 2))##
## ###############################################
## # 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.80266 -0.21963 0.01558 0.28686 0.73058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.571569 18.067986 1.913 0.0595 .
## z.lag.1 -0.037139 0.019458 -1.909 0.0601 .
## tt 0.014646 0.007209 2.032 0.0457 *
## z.diff.lag1 0.928088 0.107620 8.624 7.02e-13 ***
## z.diff.lag2 -0.251322 0.112917 -2.226 0.0290 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3849 on 76 degrees of freedom
## Multiple R-squared: 0.597, Adjusted R-squared: 0.5758
## F-statistic: 28.15 on 4 and 76 DF, p-value: 2.378e-14
##
##
## Value of test-statistic is: -1.9087 3.804 2.0874
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.47summary(ur.df(diff(Canada[, "e"]),
type = "drift", lags = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.91838 -0.24795 -0.02067 0.27069 0.74671
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.14088 0.05309 2.653 0.00965 **
## z.lag.1 -0.35887 0.07956 -4.511 2.25e-05 ***
## z.diff.lag 0.31844 0.10790 2.951 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3902 on 78 degrees of freedom
## Multiple R-squared: 0.2221, Adjusted R-squared: 0.2022
## F-statistic: 11.14 on 2 and 78 DF, p-value: 5.568e-05
##
##
## Value of test-statistic is: -4.5107 10.1751
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86summary(ur.df(Canada[, "U"],
type = "drift", lags = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72032 -0.20517 -0.04822 0.16036 1.27685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.51542 0.23809 2.165 0.0334 *
## z.lag.1 -0.05561 0.02505 -2.220 0.0293 *
## z.diff.lag 0.59746 0.09111 6.557 5.18e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3539 on 79 degrees of freedom
## Multiple R-squared: 0.3598, Adjusted R-squared: 0.3436
## F-statistic: 22.2 on 2 and 79 DF, p-value: 2.235e-08
##
##
## Value of test-statistic is: -2.2201 2.4762
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86summary(ur.df(diff(Canada[, "U"]),
type = "none", lags = 0))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.84654 -0.22553 -0.05021 0.15609 1.31737
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.43474 0.09155 -4.749 8.71e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3603 on 81 degrees of freedom
## Multiple R-squared: 0.2178, Adjusted R-squared: 0.2081
## F-statistic: 22.55 on 1 and 81 DF, p-value: 8.711e-06
##
##
## Value of test-statistic is: -4.7488
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.6 -1.95 -1.61summary(ur.df(Canada[, "rw"],
type = "trend", lags = 4))##
## ###############################################
## # 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
## -2.56910 -0.49348 -0.05407 0.53899 2.67852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.17108 9.38785 2.149 0.03503 *
## z.lag.1 -0.04740 0.02305 -2.056 0.04342 *
## tt 0.03090 0.02045 1.511 0.13521
## z.diff.lag1 0.17198 0.10812 1.591 0.11608
## z.diff.lag2 -0.02622 0.10898 -0.241 0.81056
## z.diff.lag3 -0.08185 0.10919 -0.750 0.45594
## z.diff.lag4 0.32104 0.10758 2.984 0.00388 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8257 on 72 degrees of freedom
## Multiple R-squared: 0.3609, Adjusted R-squared: 0.3076
## F-statistic: 6.776 on 6 and 72 DF, p-value: 1.009e-05
##
##
## Value of test-statistic is: -2.0558 3.2715 3.5018
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.47summary(ur.df(diff(Canada[, "rw"]),
type = "drift", lags = 3))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6181 -0.5280 -0.1168 0.5822 3.1184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2648 0.1632 1.623 0.108929
## z.lag.1 -0.3456 0.1317 -2.624 0.010540 *
## z.diff.lag1 -0.4140 0.1330 -3.112 0.002641 **
## z.diff.lag2 -0.3790 0.1210 -3.132 0.002489 **
## z.diff.lag3 -0.3979 0.1036 -3.841 0.000257 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8532 on 74 degrees of freedom
## Multiple R-squared: 0.4545, Adjusted R-squared: 0.425
## F-statistic: 15.41 on 4 and 74 DF, p-value: 3.255e-09
##
##
## Value of test-statistic is: -2.6244 3.8004
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86summary(ur.df(diff(Canada[, "rw"]),
type = "drift", lags = 0))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9840 -0.5713 -0.1145 0.5992 4.1143
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5529 0.1476 3.747 0.000336 ***
## z.lag.1 -0.5687 0.1015 -5.603 2.89e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9493 on 80 degrees of freedom
## Multiple R-squared: 0.2818, Adjusted R-squared: 0.2728
## F-statistic: 31.39 on 1 and 80 DF, p-value: 2.892e-07
##
##
## Value of test-statistic is: -5.6027 15.7327
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86#add
x1 <- summary(ur.df(Canada[, "prod"],
type = "trend", lags = 2))
x2 <- summary(ur.df(diff(Canada[, "prod"]),
type = "drift", lags = 1))
print(x1) #table 8.12 row 1 ##
## ###############################################
## # 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
## -2.19924 -0.38994 0.04294 0.41914 1.71660
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.415228 15.309403 1.987 0.0506 .
## z.lag.1 -0.075791 0.038134 -1.988 0.0505 .
## tt 0.013896 0.006422 2.164 0.0336 *
## z.diff.lag1 0.284866 0.114359 2.491 0.0149 *
## z.diff.lag2 0.080019 0.116090 0.689 0.4927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6851 on 76 degrees of freedom
## Multiple R-squared: 0.1354, Adjusted R-squared: 0.08993
## F-statistic: 2.976 on 4 and 76 DF, p-value: 0.02438
##
##
## Value of test-statistic is: -1.9875 2.3 2.3817
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.04 -3.45 -3.15
## phi2 6.50 4.88 4.16
## phi3 8.73 6.49 5.47print(x2) #table 8.12 row 2##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.05124 -0.39530 0.07819 0.41109 1.75129
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.11534 0.08029 1.437 0.155
## z.lag.1 -0.68893 0.13350 -5.160 1.83e-06 ***
## z.diff.lag -0.04274 0.11275 -0.379 0.706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6971 on 78 degrees of freedom
## Multiple R-squared: 0.3615, Adjusted R-squared: 0.3451
## F-statistic: 22.08 on 2 and 78 DF, p-value: 2.526e-08
##
##
## Value of test-statistic is: -5.1604 13.3184
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86#etc2.8.10 Canada VAR: Lag-order selection
VARselect(Canada, lag.max = 8, type = "both") #table 8-13## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 2 1 3
##
## $criteria
## 1 2 3 4 5
## AIC(n) -6.272579064 -6.636669705 -6.771176872 -6.634609210 -6.398132246
## HQ(n) -5.978429449 -6.146420347 -6.084827770 -5.752160366 -5.319583658
## SC(n) -5.536558009 -5.409967947 -5.053794411 -4.426546046 -3.699388378
## FPE(n) 0.001889842 0.001319462 0.001166019 0.001363175 0.001782055
## 6 7 8
## AIC(n) -6.307704843 -6.070727259 -6.06159685
## HQ(n) -5.033056512 -4.599979185 -4.39474903
## SC(n) -3.118280272 -2.390621985 -1.89081087
## FPE(n) 0.002044202 0.002768551 0.003060122.8.11 Diagnostic tests for VAR(p) specifications for Canadian data
Canada <- Canada[, c("prod", "e", "U", "rw")]
p1ct <- VAR(Canada, p = 1, type = "both")
p2ct <- VAR(Canada, p = 2, type = "both")
p3ct <- VAR(Canada, p = 3, type = "both")
## Serial
serial.test(p3ct, lags.pt = 16, #table 8.14 row 1 cols 2,3
type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object p3ct
## Chi-squared = 173.97, df = 208, p-value = 0.9587serial.test(p2ct, lags.pt = 16, #table 8.14 row 2 cols 2,3
type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object p2ct
## Chi-squared = 209.74, df = 224, p-value = 0.7443serial.test(p1ct, lags.pt = 16,
type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 233.5, df = 240, p-value = 0.606serial.test(p3ct, lags.pt = 16, #table 8.14 row 1 cols 4,5
type = "PT.adjusted")##
## Portmanteau Test (adjusted)
##
## data: Residuals of VAR object p3ct
## Chi-squared = 198.04, df = 208, p-value = 0.6785serial.test(p2ct, lags.pt = 16,
type = "PT.adjusted")##
## Portmanteau Test (adjusted)
##
## data: Residuals of VAR object p2ct
## Chi-squared = 236.08, df = 224, p-value = 0.2769serial.test(p1ct, lags.pt = 16,
type = "PT.adjusted")##
## Portmanteau Test (adjusted)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 256.88, df = 240, p-value = 0.2167## JB
normality.test(p3ct) #table 8.14 row 1 cols 6,7## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object p3ct
## Chi-squared = 9.665, df = 8, p-value = 0.2893
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object p3ct
## Chi-squared = 4.3714, df = 4, p-value = 0.3581
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object p3ct
## Chi-squared = 5.2936, df = 4, p-value = 0.2585normality.test(p2ct)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object p2ct
## Chi-squared = 2.2882, df = 8, p-value = 0.9709
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object p2ct
## Chi-squared = 1.2998, df = 4, p-value = 0.8614
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object p2ct
## Chi-squared = 0.9884, df = 4, p-value = 0.9115normality.test(p1ct)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 9.9189, df = 8, p-value = 0.2708
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 6.356, df = 4, p-value = 0.1741
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 3.5629, df = 4, p-value = 0.4684## ARCH
arch.test(p3ct, lags.multi = 5) #table 8.14 row 1 cols 8,9##
## ARCH (multivariate)
##
## data: Residuals of VAR object p3ct
## Chi-squared = 512.04, df = 500, p-value = 0.3451arch.test(p2ct, lags.multi = 5)##
## ARCH (multivariate)
##
## data: Residuals of VAR object p2ct
## Chi-squared = 528.14, df = 500, p-value = 0.1855arch.test(p1ct, lags.multi = 5)##
## ARCH (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 570.14, df = 500, p-value = 0.01606## Stability (Recursive CUSUM)
plot(stability(p3ct), nc = 2) #fig 8.4
plot(stability(p2ct), nc = 2) #fig 8.5
plot(stability(p1ct), nc = 2) #fig 8.6
2.8.12 Johansen cointegration tests for Canadian
summary(ca.jo(Canada, type = "trace", #table 8.15
ecdet = "trend", K = 3,
spec = "transitory"))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 4.505013e-01 1.962777e-01 1.676668e-01 4.647108e-02 2.632104e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.85 10.49 12.25 16.26
## r <= 2 | 18.72 22.76 25.32 30.45
## r <= 1 | 36.42 39.06 42.44 48.45
## r = 0 | 84.92 59.14 62.99 70.05
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.l1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000
## e.l1 -0.02385143 1.2946681 -2.8831559 4.2418087 -8.2903941
## U.l1 3.16874549 3.4036732 -7.4261514 6.8413561 -12.5578436
## rw.l1 1.83528156 -0.3330945 1.3978789 -0.1393999 2.4466500
## trend.l1 -1.30156097 -0.2302803 -0.5093218 -1.5925918 0.2831079
##
## Weights W:
## (This is the loading matrix)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.d -0.006535281 -0.02763446 -0.070975296 -0.014754352 1.077470e-11
## e.d -0.008503348 0.11414011 -0.008156659 0.003988051 7.400392e-12
## U.d -0.004718574 -0.06154306 0.020719431 -0.006557248 -4.663936e-12
## rw.d -0.046213350 -0.14579644 -0.016945105 0.011896044 6.951877e-12summary(ca.jo(Canada, type = "trace", #not shown
ecdet = "trend", K = 2,
spec = "transitory"))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 4.483918e-01 2.323995e-01 1.313250e-01 4.877895e-02 9.508809e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 4.10 10.49 12.25 16.26
## r <= 2 | 15.65 22.76 25.32 30.45
## r <= 1 | 37.33 39.06 42.44 48.45
## r = 0 | 86.12 59.14 62.99 70.05
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.l1 1.0000000 1.0000000 1.00000000 1.000000 1.000000
## e.l1 2.7132129 -6.3190324 0.49616472 16.333916 -10.368563
## U.l1 8.8369211 -15.2682881 1.48062661 25.774259 -16.048489
## rw.l1 -0.3716323 3.1817254 -0.04085215 -2.546391 4.927457
## trend.l1 -0.4177976 -0.9335588 -0.26592659 -3.413555 -1.753060
##
## Weights W:
## (This is the loading matrix)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.d 0.023155644 -0.02832697 -0.10914770 -0.006295988 -4.784843e-13
## e.d 0.005602438 -0.01739149 0.08679396 -0.001019323 -4.385463e-13
## U.d -0.019277135 0.01381763 -0.03696147 -0.002276871 4.919851e-13
## rw.d -0.084618968 -0.02739056 -0.07798404 0.003985020 -1.032563e-132.8.13 VECM with r = 1 and normalization with respect to real wages
vecm <- ca.jo(Canada[, c("rw", "prod", "e", "U")],
type = "trace", ecdet = "trend",
K = 3, spec = "transitory")
vecm.r1 <- cajorls(vecm, r = 1)
alpha <- coef(vecm.r1$rlm)[1, ]
beta <- vecm.r1$beta
resids <- resid(vecm.r1$rlm)
N <- nrow(resids)
sigma <- crossprod(resids) / N
## t-stats for alpha
alpha.se <- sqrt(solve(crossprod(
cbind(vecm@ZK %*% beta, vecm@Z1)))
[1, 1]* diag(sigma))
alpha.t <- alpha / alpha.se
## t-stats for beta
beta.se <- sqrt(diag(kronecker(solve(
crossprod(vecm@RK[, -1])),
solve(t(alpha) %*% solve(sigma)
%*% alpha))))
beta.t <- c(NA, beta[-1] / beta.se)
#add
print(beta) #table 8.16 row 1## ect1
## rw.l1 1.00000000
## prod.l1 0.54487553
## e.l1 -0.01299605
## U.l1 1.72657188
## trend.l1 -0.70918872print(beta.t) #table 8.16 row 2## [1] NA 0.90044324 -0.01917236 1.19328037 -2.56940604print(alpha) #table 8.16 row 3## rw.d prod.d e.d U.d
## -0.084814510 -0.011994081 -0.015606039 -0.008659911print(alpha.t)#table 8.16 row 4## rw.d prod.d e.d U.d
## -5.7117416 -0.9186147 -2.1579440 -1.48689892.8.14 VECM with r = 1 and normalization with respect to real wages
vecm <- ca.jo(Canada[, c("prod", "e", "U", "rw")],
type = "trace", ecdet = "trend",
K = 3, spec = "transitory")
SR <- matrix(NA, nrow = 4, ncol = 4)
SR[4, 2] <- 0
SR## [,1] [,2] [,3] [,4]
## [1,] NA NA NA NA
## [2,] NA NA NA NA
## [3,] NA NA NA NA
## [4,] NA 0 NA NALR <- matrix(NA, nrow = 4, ncol = 4)
LR[1, 2:4] <- 0
LR[2:4, 4] <- 0
LR## [,1] [,2] [,3] [,4]
## [1,] NA 0 0 0
## [2,] NA NA NA 0
## [3,] NA NA NA 0
## [4,] NA NA NA 0svec <- SVEC(vecm, LR = LR, SR = SR, r = 1,
lrtest = FALSE,
boot = TRUE, runs = 100)
svec##
## SVEC Estimation Results:
## ========================
##
##
## Estimated contemporaneous impact matrix:
## prod e U rw
## prod 0.58402 0.07434 -0.152578 0.06900
## e -0.12029 0.26144 -0.155096 0.08978
## U 0.02526 -0.26720 0.005488 0.04982
## rw 0.11170 0.00000 0.483771 0.48791
##
## Estimated long run impact matrix:
## prod e U rw
## prod 0.7910 0.0000 0.0000 0
## e 0.2024 0.5769 -0.4923 0
## U -0.1592 -0.3409 0.1408 0
## rw -0.1535 0.5961 -0.2495 0svec$SR / svec$SRse## prod e U rw
## prod 6.2282796 0.6910641 -0.6202902 1.000083
## e -1.9446274 4.4100368 -0.7922953 2.346319
## U 0.5020004 -6.5884656 0.1061060 1.619004
## rw 0.7110005 NaN 0.6947814 5.904147svec$LR / svec$LRse## prod e U rw
## prod 4.8823949 NaN NaN NaN
## e 0.9811667 3.140225 -0.7800915 NaN
## U -1.5641837 -3.781391 0.8505022 NaN
## rw -0.9278494 3.774902 -0.8493459 NaN#add
print(svec) #table 8.17##
## SVEC Estimation Results:
## ========================
##
##
## Estimated contemporaneous impact matrix:
## prod e U rw
## prod 0.58402 0.07434 -0.152578 0.06900
## e -0.12029 0.26144 -0.155096 0.08978
## U 0.02526 -0.26720 0.005488 0.04982
## rw 0.11170 0.00000 0.483771 0.48791
##
## Estimated long run impact matrix:
## prod e U rw
## prod 0.7910 0.0000 0.0000 0
## e 0.2024 0.5769 -0.4923 0
## U -0.1592 -0.3409 0.1408 0
## rw -0.1535 0.5961 -0.2495 0print(svec$LR) #table 8.18## prod e U rw
## prod 0.7910152 0.0000000 0.0000000 0
## e 0.2024150 0.5768610 -0.4922935 0
## U -0.1592277 -0.3408997 0.1408076 0
## rw -0.1534562 0.5960848 -0.2495122 02.8.15 SVEC: Overidentification test
LR[3, 3] <- 0
LR## [,1] [,2] [,3] [,4]
## [1,] NA 0 0 0
## [2,] NA NA NA 0
## [3,] NA NA 0 0
## [4,] NA NA NA 0svec.oi <- SVEC(vecm, LR = LR, SR = SR, r = 1,
lrtest = TRUE, boot = FALSE)
svec.oi <- update(svec, LR = LR, lrtest = TRUE,
boot = FALSE)
svec.oi$LRover##
## LR overidentification
##
## data: vecm
## Chi^2 = 6.0745, df = 1, p-value = 0.013712.8.16 SVEC: Impulse response analysis
svec.irf <- irf(svec, response = "U",
n.ahead = 48, boot = TRUE)
#svec.irf #commented as it's a bit long
plot(svec.irf) #fig 8.7 only U->U is switched with e->U
2.8.17 Forecast error variance decomposition of Canadian unemployment
fevd.U <- fevd(svec, n.ahead = 48)$U #in the book table 8.19 shows selected quarters (rows) from this table
#add
print(fevd.U) #table 8.19## prod e U rw
## [1,] 0.008557524 0.9577457 0.0004040634 0.033292665
## [2,] 0.003498300 0.9144309 0.0695984727 0.012472368
## [3,] 0.001602381 0.8259016 0.1668380265 0.005657955
## [4,] 0.006573055 0.7791191 0.2084044174 0.005903410
## [5,] 0.020145320 0.7474405 0.2249732969 0.007440907
## [6,] 0.033909976 0.7229202 0.2344049017 0.008764882
## [7,] 0.045028156 0.7058824 0.2388375078 0.010251975
## [8,] 0.054142944 0.6948782 0.2397585724 0.011220315
## [9,] 0.061461190 0.6880310 0.2389357068 0.011572130
## [10,] 0.067140752 0.6840169 0.2372250509 0.011617269
## [11,] 0.071601530 0.6818522 0.2350545482 0.011491762
## [12,] 0.075225208 0.6808251 0.2327031161 0.011246593
## [13,] 0.078274116 0.6804919 0.2302922101 0.010941793
## [14,] 0.080936114 0.6805744 0.2278711130 0.010618335
## [15,] 0.083332770 0.6808933 0.2254807404 0.010293217
## [16,] 0.085534118 0.6813460 0.2231450971 0.009974818
## [17,] 0.087583505 0.6818787 0.2208696581 0.009668146
## [18,] 0.089508187 0.6824605 0.2186560898 0.009375228
## [19,] 0.091323587 0.6830728 0.2165070312 0.009096551
## [20,] 0.093039536 0.6837047 0.2144237355 0.008832053
## [21,] 0.094664311 0.6843484 0.2124059956 0.008581334
## [22,] 0.096205325 0.6849976 0.2104532484 0.008343802
## [23,] 0.097669206 0.6856473 0.2085647106 0.008118798
## [24,] 0.099062088 0.6862931 0.2067392168 0.007905627
## [25,] 0.100389656 0.6869315 0.2049752782 0.007703573
## [26,] 0.101657048 0.6875598 0.2032711839 0.007511925
## [27,] 0.102868832 0.6881761 0.2016250690 0.007329994
## [28,] 0.104029047 0.6887789 0.2000349773 0.007157121
## [29,] 0.105141261 0.6893671 0.1984989084 0.006992681
## [30,] 0.106208637 0.6899404 0.1970148508 0.006836090
## [31,] 0.107233994 0.6904984 0.1955808105 0.006686804
## [32,] 0.108219868 0.6910410 0.1941948317 0.006544318
## [33,] 0.109168553 0.6915683 0.1928550086 0.006408169
## [34,] 0.110082144 0.6920804 0.1915594924 0.006277928
## [35,] 0.110962559 0.6925777 0.1903064952 0.006153202
## [36,] 0.111811570 0.6930605 0.1890942937 0.006033632
## [37,] 0.112630812 0.6935291 0.1879212289 0.005918886
## [38,] 0.113421803 0.6939838 0.1867857069 0.005808660
## [39,] 0.114185954 0.6944252 0.1856861978 0.005702677
## [40,] 0.114924580 0.6948535 0.1846212346 0.005600679
## [41,] 0.115638909 0.6952692 0.1835894121 0.005502432
## [42,] 0.116330088 0.6956728 0.1825893853 0.005407720
## [43,] 0.116999196 0.6960646 0.1816198672 0.005316342
## [44,] 0.117647240 0.6964450 0.1806796271 0.005228115
## [45,] 0.118275170 0.6968145 0.1797674887 0.005142868
## [46,] 0.118883879 0.6971733 0.1788823279 0.005060444
## [47,] 0.119474210 0.6975220 0.1780230709 0.004980698
## [48,] 0.120046955 0.6978609 0.1771886918 0.004903492