data(econ5, package = "astsa")
macro <- window(econ5, end = c(1988, 2))
x <- cbind(growth = 100 * diff(log(macro[, "gnp"])),
du = diff(macro[, "unemp"]))Workbook 13 — Fitting Time Series Models
Regression, Innovations, and Error Correction
\[\def\E{\mathbb{E}} \def\P{\mathbb{P}} \def\R{\mathbb{R}} \def\Z{\mathbb{Z}} \def\N{\mathbb{N}} \def\cH{\mathcal{H}} \def\Var{\operatorname{Var}} \def\Cov{\operatorname{Cov}} \def\Corr{\operatorname{Corr}}\]
Fitting a VAR with Ordinary Least Squares
To fit a VAR(\(p\)) to \(d\) series, we fit one multiple regression for each series. Its current value is the response; the preceding \(p\) values of all \(d\) series are the \(pd\) predictors. Each regression uses the same predictors, with its own intercept and coefficients. Collecting the equations gives the vector model \[ X_t=c+A_1X_{t-1}+\cdots+A_pX_{t-p}+Z_t, \] where \(c\) collects the intercepts, \(A_k\) collects the coefficients on lag \(k\), and \(Z_t\) collects the current errors.
Given observed vectors \(x_1,\ldots,x_T\), we pair the response \(x_{i,t}\) from series \(i\) with every coordinate of \(x_{t-1},\ldots,x_{t-p}\). Times \(t=p+1,\ldots,T\) give \(T-p\) regression observations; the first \(p\) vectors supply the initial predictors. The equation being fitted is \[ X_{i,t}=c_i+\sum_{k=1}^p\sum_{j=1}^d (A_k)_{ij}X_{j,t-k}+Z_{i,t}. \] The coefficient on series \(j\) at lag \(k\) is \((A_k)_{ij}\). Thus the coefficients for response series \(i\) form row \(i\) of each matrix \(A_k\).
The errors in different equations may be correlated at the same time. We assume i.i.d. innovation vectors with mean zero and covariance \[ \Cov(Z_t)=\Sigma_Z, \] where \(\Sigma_Z\) is an unknown, finite positive-definite matrix that is the same at every time. Its diagonal entries are innovation variances; its off-diagonal entries are covariances between errors in different equations at the same time. Those entries need not be zero. The independence assumption concerns innovation vectors at different times.
Ordinary least squares minimizes the sum of squared residuals in each response equation. Together, these regressions minimize \[ \sum_{t=p+1}^T \left\|x_t-c-\sum_{k=1}^p A_kx_{t-k}\right\|_2^2. \] The squared vector norm adds the squared residuals across the \(d\) equations. Each equation has its own coefficients, so minimizing its sum does not constrain the other equations. If the predictor columns, including the intercept, are linearly independent, the \(d\) OLS fits uniquely determine \(\widehat c\) and \(\widehat A_1,\ldots,\widehat A_p\).
Each equation has its own unrestricted coefficients and uses the same predictors. Consequently, using \(\Sigma_Z^{-1}\) to weight the residual vectors gives the same OLS coefficient estimates, even when the errors are correlated across equations. \(\Sigma_Z\) does enter the covariances between coefficient estimates from different equations. We estimate it from residual squares and cross-products after fitting.
To use these regressions for prediction, we assume that the current error has mean zero given the past: \(\E[Z_t\mid\mathcal H_{t-1}]=0\), where \(\mathcal H_{t-1}\) contains the observations through time \(t-1\). The lagged predictors are already known given that history, so the regression’s linear expression is the conditional mean of \(X_t\). After fitting, we insert the most recent \(p\) observations to estimate the next conditional mean: \[ \widehat x_{T+1\mid T} =\widehat c+\widehat A_1x_T+\cdots+\widehat A_px_{T+1-p}. \] The coefficient estimates and the observed lags together determine this one-step forecast.
Output Growth and Unemployment
Example 13.1 (Predicting unemployment changes). When production grows, firms often need more labor. The negative empirical association between output growth and changes in unemployment is known as Okun’s law. Does recent output growth also help predict the next change in unemployment, after allowing for unemployment’s own past?
The dataset astsa::econ5 contains quarterly U.S. real gross national product (GNP), in billions of 1982 dollars, and unemployment as a percentage of the labor force. We use 1948 Q3 through 1988 Q2. Write \(Q_t\) for real GNP and \(R_t\) for the unemployment rate, and model \[
X_t=\begin{pmatrix}G_t\\D_t\end{pmatrix},
\qquad G_t=100\{\log Q_t-\log Q_{t-1}\},
\qquad D_t=R_t-R_{t-1}.
\] Here \(G_t\) approximates quarterly percentage growth, without annualizing; \(D_t\) is a change in percentage points. A move in unemployment from 5% to 6% gives \(D_t=1\).
Fit the Regression Equations
We predict this quarter’s unemployment change from growth and unemployment changes in each of the preceding two quarters. These are the response and four predictors in the unemployment equation: \[ \begin{aligned} D_t={}&c_2+(A_1)_{21}G_{t-1}+(A_1)_{22}D_{t-1}\\ &\quad +(A_2)_{21}G_{t-2}+(A_2)_{22}D_{t-2}+Z_{2,t}. \end{aligned} \] Predicting growth from the same four predictors gives the other regression. Together they form a VAR(2), with the growth coefficients in the first rows of \(A_1,A_2\) and the unemployment coefficients in the second rows. The following code pairs responses with their predictors and displays the first three observations for the unemployment regression.
# embed orders columns as current vector, lag 1, then lag 2.
paired <- as.data.frame(embed(x, 3))
names(paired) <- c("growth", "du", "growth_l1", "du_l1",
"growth_l2", "du_l2")
knitr::kable(head(paired[, -1], 3), digits = 3, row.names = FALSE)| du | growth_l1 | du_l1 | growth_l2 | du_l2 |
|---|---|---|---|---|
| 0.6 | 0.821 | 0.1 | 0.918 | -0.3 |
| 1.0 | -1.171 | 0.6 | 0.821 | 0.1 |
| 1.4 | -0.586 | 1.0 | -1.171 | 0.6 |
lm fits the unemployment equation from these columns. vars::VAR constructs the lagged predictors and fits both equations by OLS:
ols_u <- lm(du ~ growth_l1 + du_l1 + growth_l2 + du_l2, data = paired)
fit_var <- vars::VAR(x, p = 2, type = "const")
coef(summary(fit_var$varresult$du)) Estimate Std. Error t value Pr(>|t|)
growth.l1 -0.21471 0.03663 -5.8613 2.755e-08
du.l1 0.02696 0.09378 0.2875 7.741e-01
growth.l2 -0.18450 0.04024 -4.5849 9.420e-06
du.l2 -0.21070 0.08148 -2.5858 1.066e-02
const 0.33536 0.05355 6.2623 3.700e-09The lm call and fit_var$varresult$du give the same coefficients. The coefficient-table row growth.l1 estimates \((A_1)_{21}\): the coefficient on the first growth lag in the unemployment equation. Its estimate is -0.215. Holding the other three predictors fixed, one percentage point more output growth in the preceding quarter lowers the fitted unemployment change by about 0.215 percentage points.
The standard error measures uncertainty in this coefficient estimate. It combines the residual variance with the spread and correlation of the lagged predictors, using the usual OLS covariance formula. The regression observations overlap in time, however: a response at one time becomes a predictor at a later time. We therefore need a time-series justification for using that formula.
For these standard errors, we use the stable causal VAR with the i.i.d. innovations specified above. Causality and independence across time make the current innovation independent of past predictors, so \[ \Cov(Z_t\mid\mathcal H_{t-1})=\Sigma_Z. \] Thus its covariance remains constant even after conditioning on the predictors. Stability ensures that past innovations have diminishing effects; with the finite moments, this allows the regression averages to converge. Together these conditions give consistent OLS estimates and justify the usual coefficient covariance formula in large samples. Gaussianity is not needed.
confint uses the estimate plus or minus a Student \(t\) critical value times its standard error, with regression residual degrees of freedom; the coefficient p-values use the same \(t\) reference. These are large-sample approximations for the time-series regression, even under Gaussian innovations.
ci_growth <- confint(fit_var$varresult$du, "growth.l1")
ci_growth 2.5 % 97.5 %
growth.l1 -0.2871 -0.1423The approximate 95% interval is [-0.287, -0.142] percentage points of unemployment change per percentage point of growth. Both endpoints are negative. This interval concerns the first growth lag; the predictive contribution of growth also includes the second lag.
Does Growth Add Predictive Information?
The question about growth as a whole compares two conditional-mean forecasts of unemployment change: one using the histories of both series, and one using only unemployment’s history. Growth Granger-causes unemployment change in mean if adding growth’s history changes that forecast with positive probability. The direction matters: predicting unemployment from growth is a different question from predicting growth from unemployment.
For this stable VAR with i.i.d. innovations having a joint density, no Granger causality from growth to unemployment change is equivalent to \[ (A_1)_{21}=(A_2)_{21}=0. \] With both coefficients zero, the unemployment conditional mean uses only its own lagged values. Conversely, the joint-density condition rules out exact relations that would make a nonzero growth contribution already determined by unemployment’s history. Thus a nonzero growth coefficient adds predictive information. The joint test assesses the two growth lags together.
causality tests these restrictions together, accounting for the estimated variances and covariance of the two growth coefficients. It scales their Wald statistic and compares it with an approximate \(F\) reference distribution under the joint null, using the conventional VAR coefficient covariance estimate.
growth_test <- vars::causality(fit_var, cause = "growth")$Granger
growth_test
Granger causality H0: growth do not Granger-cause du
data: VAR object fit_var
F-Test = 29, df1 = 2, df2 = 304, p-value = 3e-12The p-value is 3.3e-12, giving strong evidence against the two zero restrictions. In the fitted VAR(2), recent output growth adds predictive information beyond unemployment’s own history. The negative growth coefficients associate stronger recent growth with a smaller subsequent rise, or a larger fall, in unemployment. This is a predictive conclusion. An effect of a policy intervention would require assumptions about how that intervention changes the system.
The test uses the innovation assumptions of the fitted VAR. One implication of i.i.d. innovations is zero serial correlation. The following check combines residual autocorrelations and cross-correlations through lag 12 and compares its statistic with an approximate chi-squared reference distribution. Workbook 14 develops this joint test:
macro_check <- vars::serial.test(fit_var, lags.pt = 12,
type = "PT.adjusted")$serial
macro_check
Portmanteau Test (adjusted)
data: Residuals of VAR object fit_var
Chi-squared = 47, df = 40, p-value = 0.2Its p-value is 0.22, so this check finds no evidence of remaining autocorrelation through lag 12 at the 5% level.
Recover the Mean and Innovation Covariance
The const coefficients are equation intercepts. The lagged terms also contribute to the stationary mean: taking expectations in the VAR gives \(\mu=c+A_1\mu+A_2\mu\). Hence \(c=(I_2-A_1-A_2)\mu\), and we solve this equation using the fitted coefficients to estimate \(\mu\). In each coefficient matrix, rows correspond to responses and columns to lagged predictors, in the order of x:
A_hat <- vars::Acoef(fit_var)
c_hat <- vapply(fit_var$varresult, function(eq) coef(eq)["const"], numeric(1))
mu_hat <- solve(diag(2) - A_hat[[1]] - A_hat[[2]], c_hat)
names(mu_hat) <- colnames(x)
mu_hat growth du
0.80579 0.01156 vars::roots(fit_var)[1] 0.6116 0.6116 0.3624 0.3624The reported companion-matrix eigenvalue moduli are all below one, so the fitted model is stable. The entries of mu_hat are its mean quarterly growth and mean unemployment change, in the units defined above.
The remaining one-step errors estimate the innovation covariance: \[
\widehat z_t=x_t-\widehat c-\widehat A_1x_{t-1}
-\widehat A_2x_{t-2}.
\] summary(fit_var)$covres averages their squares and cross-products with a regression degrees-of-freedom adjustment. Its diagonal entries estimate the two innovation variances; the off-diagonal entry estimates their same-quarter covariance.
Sigma_hat <- summary(fit_var)$covres
Sigma_hat growth du
growth 1.0151 -0.2136
du -0.2136 0.1479This covariance describes variation left in the two response equations. The earlier interval instead used the standard error of an estimated coefficient, whose squared value is stored in vcov(fit_var$varresult$du)["growth.l1", "growth.l1"]. \(\diamond\)
Fitting a VARMA in Two Stages
A VARMA adds past innovations to the predictors used by a VAR: \[ X_t=c+\sum_{i=1}^p A_iX_{t-i}+Z_t+\sum_{j=1}^q C_jZ_{t-j}. \] For a \(d\)-dimensional series, each response equation would have \(pd+qd\) predictors and an intercept. The observation lags are available, but the innovation lags are not: the data contain \(x_t\), not \(z_t\). Computing those errors requires the coefficients we are trying to estimate. We first need a way to approximate the missing predictors.
We first fit a VAR with many observation lags and use its prediction errors to approximate the innovations. Invertibility explains why this can work. For a stationary causal, invertible VARMA with mean zero, the inverse filter from Workbook 12 gives an infinite autoregression, \[ X_t=\sum_{k=1}^{\infty}F_kX_{t-k}+Z_t. \] With i.i.d. mean-zero innovations of finite positive-definite covariance, the remaining error \(Z_t\) has conditional mean zero given past observations. Thus the sum is the conditional-mean prediction. The matrices \(F_k\) account for both AR and MA effects and need not vanish after lag \(p\): approximating this prediction can require many observation lags even when the original VARMA has only a few lags. This argument does not require Gaussian innovations.
The first stage fits a finite VAR with \(m\) observation lags. Its residuals are \[ \widetilde z_t=x_t-\sum_{k=1}^m\widehat F_kx_{t-k}. \] If these fitted predictions approximate the infinite-lag predictions, then \(\widetilde z_t\) approximates the innovation \(z_t\). The approximation has two sources of error: omitted lags beyond \(m\), and estimation of the retained coefficients. More lags reduce truncation but require estimating more coefficients. For an unknown nonzero mean, include an intercept in this preliminary regression.
The second stage treats these residuals as supplied predictors. Regress \(x_t\) on \(x_{t-1},\ldots,x_{t-p}\) and \(\widetilde z_{t-1},\ldots,\widetilde z_{t-q}\), including an intercept unless it is known to be zero. The observation-lag coefficients estimate \(A_1,\ldots,A_p\); the residual-lag coefficients estimate \(C_1,\ldots,C_q\). Both stages are ordinary least-squares calculations.
Even knowing the innovations would not always let us estimate separate AR and MA coefficients. Consider a mean-zero VARMA(1,1) with \(C=-A\): \[ X_t=AX_{t-1}+Z_t-AZ_{t-1}, \qquad X_t-Z_t=A(X_{t-1}-Z_{t-1}). \] If all eigenvalues of \(A\) lie inside the unit circle, the stationary causal solution is \(X_t=Z_t\). Thus \(A=0\), \(A=0.5I_d\), and many other choices, each paired with \(C=-A\), describe exactly the same white-noise process. These representations are all stable and invertible. The observations cannot distinguish their separate AR and MA coefficients; this is a failure of identification.
In this case, the second-stage predictor blocks \(X_{t-1}\) and \(Z_{t-1}\) are identical. The regression can estimate their combined coefficient \(A+C\), whose true value is zero, but cannot estimate \(A\) and \(C\) separately. More observations do not resolve the duplication. Replacing the innovations by estimated residuals may produce numerically distinct predictor columns, but that estimation error supplies no information about the separate matrices.
Near cancellation can also cause trouble in a finite sample: the two predictor blocks can be nearly redundant, so substantially different AR and MA estimates can give very similar fitted predictions. Accurate prediction alone therefore does not show that individual coefficients have been estimated precisely.
For the exact cancellation above, we use the white-noise model \(X_t=Z_t\) and omit the redundant AR and MA terms. More general VARMA models can have redundant combinations of matrix coefficients. A suitable set of coefficient restrictions selects one representation from those that describe the same process. Merely choosing \(p,q\) and checking stability and invertibility does not guarantee unique coefficients; see Wilms et al. (2019, Section 2).
Example 13.2 (Estimating the hidden regressors). Take \(d=2\), \(p=q=1\), and \(c=0\), writing \(A=A_1\) and \(C=C_1\). Reuse the matrices and innovation covariance from Example 12.6, taking iid Gaussian innovations for this simulation: \[ A=0.8\begin{pmatrix}\cos(\pi/6)&-\sin(\pi/6)\\ \sin(\pi/6)&\cos(\pi/6)\end{pmatrix},\qquad C=\begin{pmatrix}0.3&0.4\\0&-0.2\end{pmatrix},\qquad Z_t\overset{\mathrm{iid}}{\sim}N_2\!\left(0, \begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}\right). \]
Here \(A+C\) is lower triangular with positive diagonal entries, so it is nonsingular. This is sufficient to identify \(A\) and \(C\) within the stable, invertible VARMA(1,1) family. To see why, the first two coefficients in its infinite autoregression satisfy \[ F_1=A+C,\qquad F_2=-CF_1. \] The observation process determines these conditional-mean prediction coefficients. When \(F_1\) is nonsingular, they in turn determine \[ C=-F_2F_1^{-1},\qquad A=F_1-C. \] Thus the separate matrices are meaningful estimation targets in this example.
Both the stationary mean and the regression intercept are zero. We generate 500 observations after discarding 500 startup observations. Saving the true innovations lets us assess how well the preliminary residuals recover them.
Generate the VARMA sample (seed 54213)
set.seed(54213)
n <- 500L # retained sample size T
burn <- 500L
A <- 0.8 * matrix(c(cos(pi/6), -sin(pi/6), sin(pi/6), cos(pi/6)),
2, byrow = TRUE)
C <- matrix(c(.3, .4, 0, -.2), 2, byrow = TRUE)
Sigma <- matrix(c(1, .6, .6, 1), 2)
z <- matrix(rnorm(2 * (n + burn)), ncol = 2) %*% chol(Sigma)
y <- matrix(0, n + burn, 2)
for (tt in 2:nrow(y))
y[tt, ] <- A %*% y[tt - 1, ] + z[tt, ] + C %*% z[tt - 1, ]
y <- y[burn + seq_len(n), ]
z <- z[burn + seq_len(n), ]
colnames(y) <- c("first", "second")Take \(m=8\) for the preliminary VAR. Each equation then uses 16 observation-lag coefficients and 492 regression rows. Exercise 13.2 compares the resulting estimate with preliminary orders 4 and 12.
preliminary <- vars::VAR(y, p = 8, type = "none")
z_tilde <- matrix(NA_real_, n, 2)
z_tilde[9:n, ] <- residuals(preliminary)
The points lie close to the diagonal in both coordinates. The correlations between true innovations and preliminary residuals are 0.991 and 0.984. Here the finite VAR recovers much of the innovation variation, making its lagged residuals useful approximations to the missing predictors.
The second regression uses times 10 through 500: time 9 is the first available preliminary residual, so its lag first enters at time 10. For each response coordinate, the four predictors at time \(t\) are \(x_{1,t-1},x_{2,t-1},\widetilde z_{1,t-1},\widetilde z_{2,t-1}\). The first two coefficients form a row of \(A\), and the last two form a row of \(C\). The lm coefficient matrix puts predictors in rows and responses in columns, so the code transposes each two-row block to recover our matrix convention.
ii <- 10:n
regressors <- cbind(y[ii - 1, ], z_tilde[ii - 1, ])
stage2 <- lm(y[ii, ] ~ regressors - 1)
A_two <- t(coef(stage2)[1:2, ])
C_two <- t(coef(stage2)[3:4, ])
A_two regressorsfirst regressorssecond
first 0.6325 -0.4227
second 0.3665 0.7050C_two regressors regressors
first 0.2760 0.5176
second -0.1389 -0.1013Fitting only a VAR(1) to these observations omits \(CZ_{t-1}\). That omitted term is generally correlated with \(X_{t-1}\), which also contains \(Z_{t-1}\). The resulting OLS coefficients need not estimate \(A\).
A_omit <- vars::Acoef(vars::VAR(y, p = 1, type = "none"))[[1]]
knitr::kable(data.frame(
Fit = c("True value", "VAR(1): MA term omitted", "Two regression stages"),
A12 = c(A[1, 2], A_omit[1, 2], A_two[1, 2])), digits = 3)| Fit | A12 |
|---|---|
| True value | -0.400 |
| VAR(1): MA term omitted | -0.253 |
| Two regression stages | -0.423 |
The VAR(1) estimate of \(A_{12}\) is -0.253, compared with the true value \(-0.400\). Including the estimated innovation lag gives -0.423, much closer to the truth in this sample. The comparison illustrates the consequence of leaving the correlated MA term out of the regression.
The two stages give preliminary estimates of the VARMA coefficients without specifying an innovation density. The second-stage lm standard errors treat the estimated innovation lags as observed predictors; they do not account for the first-stage estimation. We use this calculation for coefficient estimates, rather than interpreting its ordinary regression p-values as VARMA tests. \(\diamond\)
Reading a VARMA Fit
The two-stage fit holds preliminary innovation estimates fixed during the second regression. A likelihood fit instead evaluates each proposed set of VARMA coefficients using the predictions and errors that those coefficients imply. To assign a density to the observations, we now add a distributional assumption to the stationary causal, invertible model: \[ X_t=c+\sum_{i=1}^p A_iX_{t-i}+Z_t+\sum_{j=1}^q C_jZ_{t-j}. \] Assume \[ Z_t\overset{\mathrm{iid}}{\sim}N_d(0,\Sigma_Z), \qquad \Sigma_Z\text{ positive definite}. \] The orders \(p,q\) are fixed; \(c\), the matrices \(A_i,C_j\), and \(\Sigma_Z\) are unknown. Independence is across time. Coordinates of the same innovation vector may be correlated, as described by \(\Sigma_Z\).
As in Gaussian linear regression, maximizing the likelihood favors small squared prediction errors. For a VARMA, we must reconstruct those errors from each proposed set of coefficients before we can evaluate the fit.
Evaluating a Proposed Fit
To calculate a prediction error, we need preceding observations and innovations. Some of that history is missing at the start of the sample. One approach treats the first \(s=\max(p,q)\) observed vectors as given and chooses \(q\) starting residuals, for example zeros. For each proposed set of coefficients, we then calculate the subsequent errors in time order: \[ e_t=x_t-c-\sum_{i=1}^p A_ix_{t-i}-\sum_{j=1}^q C_je_{t-j}, \qquad t=s+1,\ldots,T. \] Each error becomes an estimated innovation for later predictions. Changing the coefficients therefore changes both the predictions and the innovation history used to make them. Evaluating a new proposal requires repeating the recursion.
Evaluating the Gaussian innovation density at each error and adding its logarithm gives the initialized Gaussian log likelihood: \[ \ell_{\mathrm{init}} =-\frac{(T-s)d}{2}\log(2\pi) -\frac{T-s}{2}\log\det\Sigma_Z -\frac12\sum_{t=s+1}^T e_t^\top\Sigma_Z^{-1}e_t. \] For fixed \(\Sigma_Z\), only the last term depends on the coefficients. Maximizing this log likelihood therefore minimizes a sum of squared residual lengths, measured relative to the innovation covariance: \[ \sum_{t=s+1}^T e_t^\top\Sigma_Z^{-1}e_t =\sum_{t=s+1}^T\left\|\Sigma_Z^{-1/2}e_t\right\|_2^2. \] Here \(\|v\|_2^2\) adds the squared coordinates of \(v\), and \(\Sigma_Z^{-1/2}\) is the symmetric inverse square root of \(\Sigma_Z\). It standardizes the errors for their variances and same-time correlations. If \(\Sigma_Z=\sigma^2 I_d\), this criterion is simply \(\sigma^{-2}\sum_t\|e_t\|_2^2\), a constant multiple of the usual residual sum of squares.
For a VAR, conditioning on the first \(p\) observations gives residuals that are linear in the unknown coefficients. As we saw at the start of the workbook, the shared predictors and unrestricted coefficients let us minimize this criterion by fitting each equation with OLS, even when \(\Sigma_Z\) is not diagonal. For a VARMA, the lagged errors depend on the same coefficients we are fitting. The residuals are then nonlinear functions of the coefficients, so a single OLS fit does not minimize the criterion. Numerical optimization repeatedly changes the coefficients, reconstructs the errors, and compares the resulting fit.
We also need to estimate \(\Sigma_Z\). Its log-determinant term in the log likelihood can no longer be treated as constant. For each proposed set of coefficients, the likelihood-maximizing covariance is \[ \widehat\Sigma_Z=\frac{1}{T-s}\sum_{t=s+1}^T e_te_t^\top, \] provided this matrix is positive definite. After this substitution, maximizing the likelihood amounts to minimizing \(\log\det\widehat\Sigma_Z\) over the coefficients. Thus the fit uses the residual variances and cross-products together, rather than just adding unweighted squared errors. A two-stage fit can supply starting coefficients for this optimization.
Choosing starting residuals makes this an approximation to the Gaussian likelihood based only on the observed data. An exact calculation accounts for uncertainty about the missing initial history, including its effect on prediction-error covariances. MTS::VARMA uses conditional Gaussian fitting, with starting residuals rather than an exact treatment of that initial uncertainty.
Read the R Fit
MTS::VARMA computes its own preliminary VAR and residual regression estimates, then jointly optimizes the AR and MA coefficients using the conditional Gaussian criterion. The call below does not use our eight-lag estimates as its starting values.
The arguments p and q specify the numbers of observation and innovation lags. For the simulated example, take \(p=q=1\) and omit the intercept because the model has known mean zero:
fit_ma <- MTS::VARMA(y, p = 1, q = 1, include.mean = FALSE)For this VARMA(1,1), Phi is \(\widehat A\). The package uses a minus sign for its MA term, so \(\widehat C\) is -Theta:
A_ml <- fit_ma$Phi
C_ml <- -fit_ma$ThetaThe standard-error matrix secoef stores predictors in rows and responses in columns. Transposing its first two rows gives the standard errors for \(A\); transposing its next two gives those for \(C\). The following table compares both estimates with the simulation truth.
SE_A <- t(fit_ma$secoef[1:2, ])
SE_C <- t(fit_ma$secoef[3:4, ])
entry <- c("11", "12", "21", "22")
comparison <- data.frame(
Entry = c(paste0("A", entry), paste0("C", entry)),
Truth = c(t(A), t(C)),
Two_stage = c(t(A_two), t(C_two)),
Package = c(t(A_ml), t(C_ml)),
SE = c(t(SE_A), t(SE_C)))
knitr::kable(comparison, digits = 3)| Entry | Truth | Two_stage | Package | SE |
|---|---|---|---|---|
| A11 | 0.693 | 0.632 | 0.638 | 0.048 |
| A12 | -0.400 | -0.423 | -0.423 | 0.048 |
| A21 | 0.400 | 0.366 | 0.385 | 0.026 |
| A22 | 0.693 | 0.705 | 0.707 | 0.026 |
| C11 | 0.300 | 0.276 | 0.281 | 0.052 |
| C12 | 0.400 | 0.518 | 0.487 | 0.057 |
| C21 | 0.000 | -0.139 | -0.155 | 0.053 |
| C22 | -0.200 | -0.101 | -0.125 | 0.049 |
For \(A_{12}\), the two-stage and package estimates both round to -0.423; likelihood fitting changes little here. For \(C_{12}\), the package estimate 0.487 is closer to the true \(0.400\) than the two-stage estimate 0.518. The \(C_{21}\) estimate moves farther from its true value of zero. Likelihood fitting optimizes the Gaussian criterion; in a finite sample, this can improve some coefficient estimates while worsening others.
The coefficient \(C_{12}\) describes the contribution of the lagged innovation in coordinate 2, \(Z_{2,t-1}\), through the MA term, holding the lagged observations and the other innovation coordinate fixed. An innovation also enters \(X_{t-1}\), so its full contribution one period later includes the AR term. Substituting the equation for \(X_{t-1}\) into \(X_t=AX_{t-1}+Z_t+CZ_{t-1}\) gives the coefficient \(A+C\) on \(Z_{t-1}\). In this simulation, \(A_{12}+C_{12}=-0.4+0.4=0\), even though \(C_{12}\) is positive.
What Supports the Reported Uncertainty?
The uncertainty calculation assumes stable and invertible dynamics. For this VARMA(1,1), the corresponding checks on the fitted matrices are the AR and MA spectral radii:
c(AR = max(Mod(eigen(A_ml)$values)),
MA = max(Mod(eigen(C_ml)$values))) AR MA
0.7836 0.2005 Both radii are below one, so the fitted model is stable and invertible.
The package estimates coefficient uncertainty from the curvature of the Gaussian log likelihood near the fitted parameters, allowing the coefficients to vary jointly. A likelihood that falls slowly as a coefficient changes indicates less precise estimation. Inverting the negative log likelihood’s curvature matrix gives an estimated covariance matrix for the coefficient estimates; the square roots of its diagonal entries give the standard errors.
Identification is essential for this calculation. In the cancellation case \(C=-A\), different coefficient pairs give the same observation distribution, so the exact likelihood cannot distinguish them. Its curvature cannot provide separate standard errors for \(A\) and \(C\). Near cancellation can leave very little curvature in some directions, producing large or numerically unstable standard errors. A numerical fit returning coefficient estimates does not by itself establish that those coefficients are identifiable.
For a correctly specified, identifiable model with stable and invertible dynamics, i.i.d. Gaussian innovations, and sufficient likelihood curvature in every parameter direction, large-sample theory approximates a coefficient estimate by a normal distribution. This gives the interval \(\widehat A_{12}\pm1.96\,\operatorname{SE}(\widehat A_{12})\). For testing \(A_{12}=0\), divide the estimate by its standard error and compare with a standard normal distribution:
A_ml[1, 2] + c(-1, 1) * 1.96 * SE_A[1, 2][1] -0.5180 -0.32832 * pnorm(-abs(A_ml[1, 2] / SE_A[1, 2]))[1] 2.229e-18MTS uses this normal approximation for its printed coefficient p-values, even though it labels the statistic a t value. This tests one coefficient within the chosen model and lag orders. It is not a general Granger-causality test: in a VARMA, forecasts from observed history combine the AR and MA matrices through the inverse filter in Workbook 12.
We can minimize the same Gaussian criterion when the innovation distribution is unknown or non-Gaussian. This is called Gaussian quasi-likelihood: it uses the model’s predictions and error covariance to choose coefficients without asserting a Gaussian innovation density.
The coefficient standard errors and normal tests above remain asymptotically valid if we retain the model and curvature conditions and replace Gaussian innovations by i.i.d. mean-zero innovations with the same positive-definite covariance and finite fourth moments. This result still assumes independence across time and constant conditional covariance. Time-varying conditional variances can require different standard errors; serial dependence can also indicate missing dynamics.
The normal approximation describes the sampling variation of the coefficient estimates. The innovations themselves may have non-Gaussian tails, so these coefficient tests do not justify Gaussian probabilities for unusually large errors. Exercise 13.4 examines this distinction.
\(\diamond\)
Fitting a Cointegrated VARMA
Cointegration allows each series to be \(I(1)\) while some linear combinations of their levels are \(I(0)\). To see how these combinations enter a VARMA fit, start with a model for the levels. Let \(V_t\) collect \(d\) individually \(I(1)\) series and write \[ V_t=c+\sum_{i=1}^p A_iV_{t-i} +Z_t+\sum_{j=1}^q C_jZ_{t-j}, \qquad p\geq1,\quad q\geq0. \] Suppose there are exactly \(r\) independent cointegrating relations, where \(0<r<d\). We will rewrite the model to show how these relations can help predict changes.
Subtract \(V_{t-1}\) from both sides and collect the lagged differences: \[ \nabla V_t=c+\Pi V_{t-1} +\sum_{k=1}^{p-1}B_k\nabla V_{t-k} +Z_t+\sum_{j=1}^q C_jZ_{t-j}, \] where \[ \Pi=\sum_{i=1}^p A_i-I_d, \qquad B_k=-\sum_{i=k+1}^p A_i. \] This is an algebraic rewriting of the levels equation. The term \(\Pi V_{t-1}\) is the only one that still involves levels rather than changes or innovations.
Cointegration supplies stationary combinations of those levels. Collect \(r\) independent cointegrating vectors in the columns of a \(d\times r\) matrix \(\beta\). We fit a model in which levels enter the change equations through \(\beta^\top V_{t-1}\). Each equation uses its own linear combination of these \(r\) quantities. Collecting those coefficients in a \(d\times r\) adjustment matrix \(\alpha\) gives \[ \Pi V_{t-1}=\alpha\beta^\top V_{t-1}, \qquad \Pi=\alpha\beta^\top. \] We require \(\alpha\) to have \(r\) independent columns, so \(\operatorname{rank}(\Pi)=r\). This is the rank restriction retained when fitting the cointegrated model.
We also choose to restrict the intercept to \(c=\alpha\rho\), where \(\rho\) is an \(r\)-vector of constants. Then the intercept and level term combine as \[ c+\Pi V_{t-1}=\alpha(\beta^\top V_{t-1}+\rho). \] Define the departures from the long-run relations by \(U_t=\beta^\top V_t+\rho\). Each coordinate of \(U_t\) is \(I(0)\); the constants allow the relations to fluctuate around nonzero levels. Substitution gives the error-correction form \[ \nabla V_t=\alpha U_{t-1} +\sum_{k=1}^{p-1}B_k\nabla V_{t-k} +Z_t+\sum_{j=1}^q C_jZ_{t-j}. \] Each column of \(\alpha\) describes how the series respond to one departure. The \(d\times d\) matrices \(B_k\) multiply recent changes, and the \(d\times d\) matrices \(C_j\) multiply past innovations. The restriction on \(c\) puts the constants inside the relations, with no separate intercept in the change equations.
A relation’s scale is arbitrary: multiplying it by a nonzero constant and dividing its adjustment coefficients by that constant leaves every prediction unchanged. We therefore fix a normalization when fitting and reporting their coefficients. Normalization fixes how the relations are written without changing the model’s predictions.
Fitting the Relations and Dynamics
If \(\beta\) and \(\rho\) were known, we could calculate the departures from observed levels. The lagged changes are observed too. For \(q>0\), however, the innovation lags are still missing. Regressing changes only on departures and lagged changes leaves out the MA contribution, which can be correlated with those predictors. Knowing the long-run relations therefore does not by itself make this an ordinary OLS fit.
For Gaussian likelihood fitting, assume \[ Z_t\overset{\mathrm{iid}}{\sim}N_d(0,\Sigma_Z), \qquad \Sigma_Z\text{ positive definite}, \] with each innovation independent of preceding observations. The unknown \(\Sigma_Z\) is constant over time and may have nonzero off-diagonal entries. We use an invertible MA representation, so past innovations can be recovered from the full observed history. The conditional-mean prediction includes all three terms: the level departures, lagged changes, and lagged innovations. The current \(Z_t\) has conditional mean zero.
In a finite sample, reconstruct the innovations for each proposed set of coefficients, as in the preceding VARMA fit. Treat the first \(s=\max(p,q)\) observations as given and choose starting residuals. Then calculate \[ \begin{aligned} e_t={}&\nabla v_t-\alpha(\beta^\top v_{t-1}+\rho)\\ &-\sum_{k=1}^{p-1}B_k\nabla v_{t-k} -\sum_{j=1}^q C_je_{t-j}, \qquad t=s+1,\ldots,T. \end{aligned} \] Changing the relations, adjustment coefficients, or short-run dynamics changes this entire residual sequence. The starting residuals approximate unobserved initial innovations; invertibility makes their effect on later errors decay. Up to a parameter-independent constant, the initialized Gaussian log likelihood is \[ \ell_{\mathrm{init}} =-\frac{T-s}{2}\log\det\Sigma_Z -\frac12\sum_{t=s+1}^T\left\|\Sigma_Z^{-1/2}e_t\right\|_2^2. \] Thus the squared-error connection remains, but the errors now depend on the unknown long-run relations as well as the AR and MA dynamics.
With the orders \(p,q\) and rank \(r\) fixed, numerical optimization changes the coefficients, reconstructs the errors, and maximizes this criterion. The fit must preserve the rank and constant restriction. It also uses an identifiable VARMA representation, as discussed above.
The Cointegrated VAR Case
When \(q=0\), the innovation lags disappear. This special case is a cointegrated VAR, usually fitted in its vector error-correction model (VECM) form. With \(\beta,\rho\) supplied, each change equation is now an ordinary regression: its response is one coordinate of \(\nabla v_t\), and its predictors are the \(r\) departures in \(u_{t-1}\) and all coordinates of \(\nabla v_{t-1},\ldots,\nabla v_{t-p+1}\). Times \(p+1,\ldots,T\) supply \(T-p\) observations. OLS estimates \(\alpha\) and the \(B_k\).
When the relations are unknown, an unrestricted regression on \(v_{t-1}\) and lagged changes generally estimates a full-rank level matrix. It need not preserve \(\Pi=\alpha\beta^\top\) with rank \(r\), or give an intercept of the required form \(\alpha\rho\). The joint fit must estimate the relations while retaining these restrictions.
For any proposed \(\beta,\rho\), calculate the departures and fit all the change equations by OLS. These regressions give the best \(\alpha,B_1,\ldots,B_{p-1}\) for that proposed relation. As before, common predictors and unrestricted coefficients in these regressions make OLS appropriate even with a nondiagonal \(\Sigma_Z\). Let \(e_t\) be their vector of residuals. The Gaussian log likelihood developed above uses their squared lengths \(e_t^\top\Sigma_Z^{-1}e_t\); estimating the covariance gives \[ \widehat\Sigma_Z=\frac{1}{T-p}\sum_{t=p+1}^T e_te_t^\top. \] The joint fit chooses the relations that minimize \(\log\det\widehat\Sigma_Z\), with the change equations refitted for each proposed set of relations. Thus the relations are chosen by how well they help predict changes in the whole system.
For \(q=0\), the Johansen procedure solves this problem by reduced-rank regression. In urca::ca.jo, the calculation first regresses both the current changes and the lagged levels on the lagged changes and keeps their residuals. The levels also include a constant column when the relations have constants. This removes the part explained by recent changes. An eigenvalue calculation then finds the \(r\) combinations of the remaining level variables most strongly associated with the remaining changes. These combinations solve the Gaussian likelihood problem with rank \(r\).
After we specify the rank, urca::cajorls takes the corresponding relations from ca.jo, normalizes their coefficients, and fits the change equations by OLS using their lagged departures. This last OLS step completes the joint fit: it uses the relations already selected by the system likelihood.
Treasury Rates and Adjustment
We illustrate the \(q=0\) case with Treasury rates, using ca.jo and cajorls. These functions fit a cointegrated VAR; they do not estimate the MA coefficients in the general model above.
A longer-term Treasury yield reflects expected future short rates and a term premium for holding the longer bond (Federal Reserve Bank of New York). This motivates looking for co-movement between short and long yields. A relation between their levels may help predict their changes, even when a model of recent changes alone misses that information.
Example 13.3 (Predicting rate changes from a level departure). The dataset urca::Raotbl1 contains natural logarithms of U.S. three-month and ten-year Treasury rates, quoted in percent per year. Keep the 143 jointly observed quarters, 1953 Q2–1988 Q4. Let \(S_t\) and \(L_t\) denote these log rates and write \(V_t=(S_t,L_t)^\top\). The supplied columns are already logged.
data(Raotbl1, package = "urca")
observed <- complete.cases(Raotbl1[, c("r3m", "r10y")])
rates <- ts(as.matrix(Raotbl1[observed, c("r3m", "r10y")]),
start = c(1953, 2), frequency = 4)
colnames(rates) <- c("short", "long")
To measure how far the short log rate is from a proposed relation \(S_t=a+bL_t\), subtract the level implied by the long log rate: \[ U_t=S_t-a-bL_t. \] For this to be a cointegrating relation under the course definition, \(S_t\) and \(L_t\) must each be \(I(1)\) and \(U_t\) must be \(I(0)\). The time-series plots alone cannot determine these properties. As an initial check, apply an ADF test to each log rate, using an intercept and an AIC-selected lag order up to four. The unit-root null is rejected when the statistic is more negative than its critical value.
Check the individual log rates
adf_checks <- sapply(1:2, function(j) {
test <- urca::ur.df(rates[, j], type = "drift", lags = 4, selectlags = "AIC")
c(Statistic = test@teststat[1], Critical_5pct = test@cval[1, "5pct"])
})
colnames(adf_checks) <- c("Short log rate", "Long log rate")
adf_checks # Reject the unit-root null when the statistic is more negative. Short log rate Long log rate
Statistic -3.171 -1.966
Critical_5pct -2.880 -2.880The short-rate statistic is -3.17, below the 5% critical value -2.88; the long-rate test does not reject. Thus these tests do not support treating both rates as \(I(1)\). We nevertheless fit a model in which this single level departure helps predict both rate changes. We impose that shared structure and estimate its coefficients, including \(a\) and \(b\). The calculation shows how to fit and read such a model; the evidence above leaves its cointegration interpretation in doubt.
For each rate, predict its current change using the preceding level departure and the two preceding changes in both rates. A positive departure means that the short log rate is above \(a+bL_t\). The two change equations together are \[ \nabla V_t= \begin{pmatrix}\alpha_S\\\alpha_L\end{pmatrix}U_{t-1} +B_1\nabla V_{t-1}+B_2\nabla V_{t-2}+Z_t. \] This is the rank-one, \(q=0\) model above, with \(\beta=(1,-b)^\top\), \(\rho=-a\), and \(p=3\). The same departure enters both equations, with adjustment coefficients \(\alpha_S\) and \(\alpha_L\). A negative \(\alpha_S\), for example, makes a positive departure lower the predicted short-rate change.
Fit the Relation, Then the Change Equation
There are two sets of coefficients to estimate: \(a,b\) determine the long-run relation, while \(\alpha_S,\alpha_L,B_1,B_2\) describe the subsequent changes. A direct two-step approach first regresses \(S_t\) on an intercept and \(L_t\). OLS chooses \(\widehat a,\widehat b\) to minimize \(\sum_t(S_t-a-bL_t)^2\). From the observations, this gives the estimated departures \[ \widehat u_t=s_t-\widehat a-\widehat b\,l_t. \] These are the residuals from the levels regression.
Now hold the estimated departures fixed. To fit the short-rate equation, regress the current change \(\nabla s_t\) on the five predictors \[ \widehat u_{t-1},\quad \nabla s_{t-1},\quad\nabla l_{t-1},\quad \nabla s_{t-2},\quad\nabla l_{t-2}. \] The first coefficient estimates \(\alpha_S\): it measures how the short rate adjusts to the preceding departure after allowing for recent changes. The other four coefficients estimate the first rows of \(B_1\) and \(B_2\). Using \(\nabla l_t\) as the response with the same predictors fits the long-rate equation.
The constant is already inside \(\widehat u_{t-1}\), so the change regression uses -1 to omit a separate intercept. The code below fits the levels relation and then the short-rate change equation.
long_run <- lm(short ~ long, data = as.data.frame(rates))
u_hat <- residuals(long_run)
dr <- diff(rates)
ii <- 3:nrow(dr)
ecm_data <- data.frame(
ds = dr[ii, "short"], dl = dr[ii, "long"],
u_l1 = u_hat[ii],
ds_l1 = dr[ii - 1, "short"], dl_l1 = dr[ii - 1, "long"],
ds_l2 = dr[ii - 2, "short"], dl_l2 = dr[ii - 2, "long"])
ecm_ols <- lm(ds ~ u_l1 + ds_l1 + dl_l1 + ds_l2 + dl_l2 - 1,
data = ecm_data)
coef(long_run)(Intercept) long
-0.6428 1.2156 coef(ecm_ols) u_l1 ds_l1 dl_l1 ds_l2 dl_l2
-0.1598 0.2482 0.8837 -0.1390 -0.3729 Because dr[ii, ] is the change from level row ii to row ii + 1, u_hat[ii] is the preceding level departure. Its fitted coefficient is -0.160. Its negative sign means that a positive departure lowers the fitted next short-rate change, holding recent changes fixed.
Read the System Fit
The levels regression chose \(a,b\) to make \(s_t-a-b l_t\) small. We now use the joint Gaussian fit to choose \(a,b\) by how well the relation helps predict both rate changes. Its estimates need not equal the levels-regression estimates: the two methods minimize different criteria.
In ca.jo, ecdet = "const" puts the constant inside the relation, and spec = "transitory" puts the level departure at time \(t-1\), as in our equation. Since \(\nabla V_{t-2}=V_{t-2}-V_{t-3}\), rewriting the model in levels uses lags through three. The argument K counts these levels lags, so we set K = 3 to include two lagged differences.
The same call computes trace tests for the rank of the lagged-level matrix. Each compares the best Gaussian likelihood under a rank restriction with a less restricted fit; large statistics argue against the restriction. The two rows below test rank zero and rank at most one. Their critical values come from the cointegration model’s large-sample theory and depend on the specification of constants and trends.
jo <- urca::ca.jo(rates, type = "trace", ecdet = "const",
K = 3, spec = "transitory")
knitr::kable(data.frame(
Null = c("Rank zero", "Rank at most one"),
Statistic = rev(jo@teststat), Critical_5pct = rev(jo@cval[, "5pct"])),
digits = 2)| Null | Statistic | Critical_5pct | |
|---|---|---|---|
| r = 0 | | Rank zero | 22.29 | 19.96 |
| r <= 1 | | Rank at most one | 4.78 | 9.24 |
At 5%, the trace tests reject rank zero and do not reject rank at most one. These results favor a rank-one level matrix for the specified lags and constant. They do not resolve the earlier short-rate result: a stationary short rate alongside an \(I(1)\) long rate can itself produce rank one, without a cointegrating relation between two \(I(1)\) rates. We retain the imposed rank for the calculation below.
cajorls with rank one uses the first relation estimated by ca.jo and normalizes the short-rate coefficient to one. It then fits both change regressions using this relation and the four lagged-change predictors. It returns the normalized relation in beta and these OLS fits in rlm.
fit_ec <- urca::cajorls(jo, r = 1)
fit_ec$beta ect1
short.l1 1.0000
long.l1 -1.2027
constant 0.6547coef(fit_ec$rlm) short.d long.d
ect1 -0.1516 0.01998
short.dl1 0.2429 -0.03678
long.dl1 0.8974 0.39218
short.dl2 -0.1456 -0.02086
long.dl2 -0.3595 -0.11029The normalization fixes the coefficient of \(S_t\) at one. The fitted relation is \[
\widehat u_t=s_t
-1.203\,l_t
+0.655.
\] In coef(fit_ec$rlm), each column is a response equation. The row ect1 contains \(\widehat\alpha_S\) and \(\widehat\alpha_L\); short.dl1 multiplies \(\nabla S_{t-1}\). The remaining lag rows are read in the same way. The constant in the displayed relation is \(-\widehat a\).
The table gives approximate 95% intervals for the adjustment coefficients under the specified cointegrated model and its i.i.d. Gaussian innovations. summary(fit_ec$rlm) computes ordinary regression standard errors with the fitted level relation held fixed. In a correctly specified cointegrated system, estimating that relation does not change the leading large-sample distribution of the adjustment estimates. This justifies using those regression standard errors here; it is a result specific to cointegration, rather than a general rule for estimated predictors.
Approximate normality gives the intervals \(\widehat\alpha\pm1.96\,\operatorname{SE}(\widehat\alpha)\). Gaussianity is sufficient here; the same coefficient approximation also holds with i.i.d. mean-zero innovations of constant positive-definite covariance and finite fourth moments. The calculation treats the rank, lag orders, and placement of the constant as fixed. It describes coefficient uncertainty under the assumed model; it does not resolve the conflicting evidence about whether both rates are \(I(1)\).
alpha_hat <- coef(fit_ec$rlm)["ect1", ]
ec_summaries <- summary(fit_ec$rlm)
alpha_se <- vapply(ec_summaries,
function(s) coef(s)["ect1", "Std. Error"], numeric(1))
knitr::kable(data.frame(Equation = c("Short log rate", "Long log rate"),
Adjustment = alpha_hat, SE = alpha_se,
Lower = alpha_hat - 1.96 * alpha_se,
Upper = alpha_hat + 1.96 * alpha_se),
row.names = FALSE, digits = 3)| Equation | Adjustment | SE | Lower | Upper |
|---|---|---|---|---|
| Short log rate | -0.152 | 0.062 | -0.274 | -0.030 |
| Long log rate | 0.020 | 0.028 | -0.035 | 0.075 |
A departure \(\widehat u_{t-1}=0.01\) means the short rate is about 1% above the level implied by the fitted relation. Holding lagged changes fixed, the adjustment term contributes approximately -0.0015 to the next short log-rate change: a fall of about 0.152 percent in the short rate. In this fitted model, most adjustment occurs through the short rate. The long-rate coefficient is small relative to its standard error, and its interval includes zero.
vec2var converts the fitted VECM to its equivalent levels VAR for methods such as residual checking. It does not fit a new model.
fit_levels <- vars::vec2var(jo, r = 1)
rates_check <- vars::serial.test(fit_levels, lags.pt = 12,
type = "PT.adjusted")$serial
rates_check
Portmanteau Test (adjusted)
data: Residuals of VAR object fit_levels
Chi-squared = 45, df = 38, p-value = 0.2The joint residual check has p-value 0.20, so it finds no evidence of remaining autocorrelation through lag 12 at 5%. The fitted rank-one model predicts that departures are followed mainly by short-rate adjustment. Deciding whether to interpret this as adjustment between cointegrated rates still requires resolving the conflicting integration evidence; a residual autocorrelation test does not address that question. \(\diamond\)
Exercises
Exercise 13.1 (Read the unemployment equation). Use Example 13.1.
Identify the matrix entry estimated by
growth.l2. Interpret its estimate and a 95% interval in the units of the two series.Is
constthe estimated mean unemployment change? Use the fitted matrix equation to obtain that mean.The p-value for
du.l1is large. Does that settle whether unemployment changes Granger-cause output growth? Name the restrictions for that question.
Exercise 13.2 (Read the innovation coefficients). Use Example 13.2.
Explain why the preliminary residuals are needed. Which regressors would be available if the true innovations were supplied?
Read \(C_{12}\) and its standard error from
fit_ma. Interpret its contribution through the MA term, respecting the package’s sign convention. Why does \(C_{12}\) alone not give the full one-period contribution of the innovation in coordinate 2 to coordinate 1?Repeat both regression stages, using 4 and then 12 lags in the preliminary VAR. Compare the resulting second-stage \(A_{12}\) estimates with the eight-lag fit. Does closer agreement with the truth in this one simulation identify a generally best preliminary order?
Exercise 13.3 (Read the adjustment). Use Example 13.3.
Translate the three entries of
fit_ec$betainto \(U_t=S_t-a-bL_t\). Which entry gives \(-a\)?Set \(u_{t-1}=0.02\). Calculate the contribution of the adjustment term to each next log-rate change, holding lagged changes fixed. Explain the signs and distinguish a percent change in a yield from a percentage-point change in a yield.
Why does a VAR fitted only to
diff(rates)leave out part of this model? What doesK = 3mean in the Johansen call?
Exercise 13.4 (A good mean fit need not give the right tails). Consider the stationary AR(1) \[
X_t=0.6X_{t-1}+Z_t,\qquad
Z_t\overset{\mathrm{iid}}{\sim}0.95N(0,1)+0.05N(0,25).
\] Generate the sample below and fit arima(x_mix, order = c(1, 0, 0), method = "ML"). Here order = c(1, 0, 0) specifies one AR lag, no differencing, and no MA lags. For this undifferenced model, arima reports the fitted mean \(\mu\) as intercept, so its equation is \(X_t-\mu=\phi(X_{t-1}-\mu)+Z_t\). The returned sigma2 estimates the innovation variance. The option method = "ML" uses Gaussian densities, so for these mixture innovations it is Gaussian quasi-likelihood. For this stationary AR(1), it consistently estimates the AR coefficient, mean, and innovation variance despite the non-Gaussian mixture density.
set.seed(542134)
z_mix <- rnorm(800, sd = ifelse(runif(800) < .05, 5, 1))
x_mix <- as.numeric(filter(z_mix, filter = .6, method = "recursive"))[301:800]Report the fitted AR coefficient, mean, and innovation variance. What are their values in the generating model?
Does an AR estimate near \(0.6\) establish that the innovations are Gaussian? Compute \(\P(|Z_t|>3)\) under the supplied mixture and under a Gaussian distribution with the same variance.
Further references: vars documentation, urca documentation, and MTS documentation; Tsay (2014), Multivariate Time Series Analysis with R and Financial Applications, Chapters 2–4, for regression, VARMA estimation and error correction. The code uses vars 1.6-1, MTS 1.2.1, urca 1.3-4 and astsa 2.5.