So here's the framework: I'm dealing with a time series model of the type \begin{align*} y_t &= a_{11} + a_{12} x_{t-1} + a_{13} z_{t-1} + a_{14} y_{t-1} + u_{1,t}\\ x_t &= a_{21} + a_{22} x_{t-1} + a_{23} z_{t-1} + u_{2,t} \\ z_t &= a_{31} + a_{33} z_{t-1} + u_{3,t} \end{align*} where
- $(x_t)$ and $(y_t)$ are cointegrated.
- $(z_t)$ is $I(0)$.
- The unobservable error term $u_{t}$ follows a stationary VAR(1) process.
That is, I'm dealing with a (restricted) VAR(1) process where some of the components are non-stationary and where the error term displays serial correlation.
I'm really only interested in estimating the parameters $(a_{11}, a_{12}, a_{13}, a_{14})$ and forecasting the $(y_t)$ series, so my question is: what is the best approach for achieving this?
Below is an R code which simulates a process of the type I'm considering. Interestingly, if I estimate the parameters via least squares it just "works", although in the literature this seems to be regarded as bad practice.
a = matrix(0, nrow = 3, ncol = 4) # VAR coefficient matrix
a[1,1] = .5
a[1,2] = .6
a[1,3] = 1
a[1,4] = .4
a[2,1] = 1
a[2,2] = 1
a[3,3] = .7
n = 200 # sample size
# generating z
sigma_u3 = 1.5
u3 = rnorm(n, mean = 0, sd = sigma_u3)
z = numeric(n)
z[1] = rnorm(1, mean = 0, sd = sigma_u3^2 / (1 - a[3,3]^2))
for (t in 2:n) {
z[t] = a[3,4]*z[t-1] + u3[t]
}
# generating x
sigma_u2 = 1
u2 = rnorm(n, mean = 0, sd = sigma_u2)
x = numeric(n)
for (t in 2:n) {
x[t] = a[2,1] + a[2,2]*x[t-1] + u2[t]
}
# generating y
sigma_e = .4
e = rnorm(n, mean = 0, sd = sigma_e)
gamma = .5 # AR(1) coefficient of the error term u1[t]
u1 = numeric(n)
u1[1] = rnorm(1, mean = 0, sd = sigma_e^2 / (1 - gamma^2))
for ( t in 2:n) {
u1[t] = gamma*u1[t-1] + e[t] # the error term u1[t] is an AR(1)
}
y = numeric(n)
for (t in 2:n) {
y[t] = a[1,1] + a[1,2]*x[t-1] + a[1,3]*z[t-1] + a[1,4]*y[t-1] + u1[t]
}
dev.new()
plot(1:n, x, type = 'l', ann = FALSE)
lines(y, col = 'red')
lines(z, col = 'blue')