Non-stationary VAR estimation 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')

 A: I suppose by "nonstationary" you mean "integrated". If so, then


*

*equation 3 must be misspecified as the the right hand side is integrated while the left hand side is stationary, so the two diverge; 

*if $x_t$ and $y_t$ are not cointegrated, equation 1 will be misspecified as its left hand side (driven by one integrated process) will be diverging from its right hand side (driven by two different integrated processes without a common stochastic trend).


Edit: After the edit of the OP, here are some news observations:


*

*$a_{12}$ and $a_{14}$ should be fixed at values yielding cointegration between $x_t$ and $a_{12} x_{t-1} + a_{14} y_{t-1}$, otherwise the left hand side will diverge from the right hand side. You can obtain these values from estimating the cointegrating vector between $x_t$ and $y_t$.

*Equation 2 will be misspecified unless $a_{22}=1$; if $a_{22} \neq 1$, the left hand side will diverge from the right hand side.


Once you fix the above two bullet points, you can estimate the equations using feasible GLS or equation-by-equation OLS. This way you will get the estimates of $a_{11}$ through $a_{41}$.
