Joint Mean-Vector and Var-Matrix estimation
for Locally Stationary VAR(1) processes

In this paper we consider $r$-dimensional multivariate data $\bm{X}_{T}(1),\dots\bm{X}_{T}(t)$ generated by a locally stationary process, and our goal is to fit to the data a parametric model with time-varying coefficients. The notation $\bm{X}_{T}(t)$ emphasizes that the data is a triangular array where at each $t$, the structure of the process depends on the sample size $T$.

To introduce the problem, consider the following uni-variate zero-mean autoregressive model of order $p$,

$$\sum_{j=0}^{p}a_{j}(\tfrac{t}{T})X_{T}(t-j)=\sigma(\tfrac{t}{T})\varepsilon(t),$$

(1)

or AR($p$), where the coefficients $a_{j}(u)$ are differentiable for $u\in(0,1)$ with bounded derivatives.

In terms of modeling, local stationarity means that if the parameters $a_{j}(u)$ are smooth in rescaled time $u$ and $T$ is large, then $a_{j}(\tfrac{s}{T})\approx a_{j}(\tfrac{t}{T})$ for all $s$ in a neightbour of $t$. For estimation, rescaling in time allows to apply non-parametric methods to recover the unknown curves. In the frequency domain, the importance of rescaling time $t$ by the sample size $T$ and developing the analysis in rescaled time $u\in(0,1)$ relies upon the uniqueness of the transfer function. Dahlhaus (1996) introduced the spectral representation of a locally stationary process

$$X_{T}(u)=\int_{-\pi}^{-\pi}\exp(i\lambda t)A_{T}^{0}(t,\lambda)d\xi(\lambda),$$

where $\xi(\lambda)$ is a stochastic process with orthogonal increments, and where the sequence $A_{T}^{0}(t,\lambda)$ converges (uniformly in $t$ and $\lambda$, as $T$ diverges) to another function $A(u,\lambda)$:

$$\sup_{t,\lambda}|A_{T}^{0}(t,\lambda)-A(\tfrac{t}{T},\lambda)|=\mathcal{O}(\tfrac{1}{T}).$$

If $A$ is smooth in $u$, then the time-varying spectral density $f(u,\lambda)=|A(\tfrac{t}{T},\lambda)|^{2}$ is uniquely determined from the triangular array.

The dichotomy between $A_{T}^{0}$ and $A$ is particularly relevant in the case of AR processes. To see this, consider the simple case where $p=1$ and and $\sigma(u)\equiv 1$. In the stationary case where the coefficient $a$ is time-invariant, the process in (1) can be written as

$$X(t)=\sum_{k=0}^{\infty}\psi_{\,\!{}_{k}}\varepsilon_{t-k}$$

with $\psi_{\,\!{}_{k}}=a^{k}$. By contrast, the locally stationary process $X_{T}(t)=a(\tfrac{t}{T})X_{T}(t-1)+\varepsilon(t)$ does not have a solution of the form

$$X_{T}(t)=\sum_{k=0}^{\infty}\psi_{\,\!{}_{k}}(u)\varepsilon_{t-k},$$

but only of the form $X_{T}(t)=\displaystyle\sum_{k=0}^{\infty}\psi_{\!{}_{T,k}}(t)\varepsilon_{t-k}$,   with   $\sup_{t}\displaystyle\sum_{k\in\mathbb{Z}}|\psi_{\!{}_{T,k}}(t)-\psi_{\,\!{}_{k}}(\tfrac{t}{T})|=\mathcal{O}(\tfrac{1}{T})$.

The seminal papers on local stationarity (Dahlhaus, 1996, 1997) provide details on the mathematics in the frequency domain. For an overview on multivariate locally stationary processes, see Dahlhaus (2012, Section 7.2).

Without loss of generality, assume that $\sigma(u)$ is known and time-invariant, that is, $\sigma(u)\equiv\sigma$. Suppose that the vector of interest $\bm{a}(u)=[a_{1}(u),\dots,a_{p}(u)]^{\top}$ depends on a finite dimensional parameter. For example, if the coefficients are polynomials in time

$\displaystyle\begin{split}a_{j}(u)=&\sum_{k=1}^{K}\theta_{jk}f_{k}(u),\quad 1\leq k\leq K,&\quad 1\leq j\leq p,\\ \mbox{with}\quad f_{k}(u)=&\,\,\,u^{k-1},\qquad\qquad 1\leq k\leq K,&\end{split}$

(2)

estimating the time-varying vector $\bm{a}(u)$ at $u\in(0,1)$ translates into estimating the time-invariant vector $\bm{\theta}=[\bm{\theta}_{1}^{\top},\dots,\bm{\theta}_{p}^{\top}]^{\top}$, where $\bm{\theta}_{j}=[\theta_{j1},\dots,\theta_{jK}]^{\top}$, with $1\leq j\leq p$.

The specification in (2) approximates the coefficients vector $\bm{a}(u)$ by global polynomials. Dahlhaus (1997, Section 4) obtained an explicit formula for the vector ${\bm{\theta}}$ as the solution of a linear system similar to the Yule-Walker equations. In the univariate setting, Dahlhaus (1996) estimates the time-varying parameters $\bm{a}(u)$ by kernel smoothers, that is, using a local-constant approximation. In the multivariate settings, Dahlhaus (2000, p. 1776) mentions the possibility of estimating the unknown parameters $\bm{\theta}(u)$ by minimazing of a local-polynomial approximation of the local likelihood.

Zhou and Wu (2010) consider univariate linear models of the form $\bm{Y}(t)=\bm{X}^{\top}\bm{\beta}(t)+\bm{\varepsilon}(t)$, where both $\bm{X}$ and $\bm{\varepsilon}$ are asumed to be locally stationary, and estimate their time-varying vector of coefficients $\bm{\beta}(t)$ by means of local polynomials.

The contribution of this paper is threefold. In terms of modeling, we consider a multivariate version of model (1) and estimate the time varying parameters in time domain. Our main contribution is the closed-form definition of the local-linear estimator of the parameters. Finally, we emphasize that the estimation of time-varying mean and time-varying AR matrix is performed jointly.

In Section 2 we derive the localized Yule-Walker equations for a locally-stationary zero-mean VAR process. In Section 3 we consider a locally-stationary VAR with time-varying mean. First, we derive the local-constant weighted-least-squares estimator, see Proposition 1. Then in Theorem 1 we establish our main result, the closed form definition of the local-linear weighted-least-squares estimator. In Section 4 we illustrate and compare the performance of the two weighted-least-squares estimators (local-constant and local-linear). Section 5 concludes and highlights the extension of our approach to the high-dimensional ($r>T$) setting.

Through the paper we use bold uppercase letters to denote matrices, and bold slanted to denote vectors. We denote by $\mathbf{I}_{m}$ the identity matrix of size $m$, by $\bm{1}_{m}$ the $m\times 1$ vector of ones, by ${\rm tr}\{\mathbf{A}\}$ the trace of $\mathbf{A}$, by $\mathbf{A}^{\top}$ the transpose of $\mathbf{A}$, by $\|\mathbf{A}\|$ the Frobenius norm $\|\mathbf{A}\|=[{\rm tr}\{\mathbf{A}^{\top}\mathbf{A}\}]^{\scalebox{0.6}{$1/2$}}$, and by $\mathbf{A}^{-1}$ the inverse of $\mathbf{A}$, that is, the square matrix such that $\mathbf{A}^{-1}\mathbf{A}=\mathbf{A}\mathbf{A}^{-1}=\mathbf{I}$. Finally, we use the acronyms VAR and WLS for vector auto-regression and weighted least squares, respectively.

Consider the following $r$-dimensional Locally Stationary Vector Auto Regression of order 1

$$\underset{r\times 1}{\bm{X}_{t}}=\underset{r\times r}{\mathbf{A}(\tfrac{t}{T})}\,\underset{r\times 1}{\bm{X}_{t-1}}+\underset{r\times 1}{\bm{\varepsilon}_{t}},\qquad t=1,\dots,T,$$

(3)

with $\bm{X}_{0}=\bm{0}$ and $\mathbb{E}[\bm{\varepsilon}_{t}\bm{X}_{s}^{\top}]=\mathbbm{1}_{\{s=t\}}\mathbf{\Gamma}_{\varepsilon}$. If the largest eigenvalue of $\mathbf{A}$ lies inside the unit circle

$$\sup_{u\in(0,1)}\left|{\rm v}_{1}[\mathbf{A}(u)]\,\right|<1,$$

$\bm{X}_{t}$ is locally stationary and causal. Our goal is to estimate $\mathbf{A}(u)$ at a fixed $u\in(0,1)$ using a localized version of the Yule-Walker equations. If we assume that the matrix-valued function $\mathbf{A}(x)$ is smooth in $x$, we can write the Taylor expansion of $\mathbf{A}(x)$ around $u$:

$$\mathbf{A}(x)=\sum_{j=0}^{\infty}\tfrac{\big{(}x-u\big{)}^{j}}{j!}\mathbf{A}^{(j)}(u)\\ =\mathbf{A}(u)+(x-u)\mathbf{A}^{(1)}(u)+\mathcal{O}([x-u]^{2})$$

where $\mathbf{A}^{(j)}(u):=\tfrac{d^{j}\mathbf{A}(x)}{dx^{j}}|_{x=u}$. We are interested in evaluating the function $\mathbf{A}(\tfrac{t}{T})$ at those value of $t$ in a neighborhood of $u$. For example, for a fixed $u_{0}\in(0,1)$ let $t_{0}=\lfloor u_{0}T\rfloor$, where $\lfloor x\rfloor$ is the largest integer not exceeding $x$. Then we have the following uniform bound:

$$\sup_{u_{0}\in(0,1)}|\tfrac{t_{0}}{T}-u_{0}|<\tfrac{1}{T}.$$

As a consequence, assuming that

$$\sup_{u\in(0,1)}\|\mathbf{A}^{(1)}(u)\|<\infty$$

we obtain following bound

$$\|\mathbf{A}(\tfrac{t}{T})-\mathbf{A}(u)\|\leq|\tfrac{t}{T}-u|\,\|\mathbf{A}^{(1)}(u)\|=\mathcal{O}(\tfrac{1}{T})\times\mathcal{O}(1)=\mathcal{O}(\tfrac{1}{T})$$

uniformly in $u$. Let $K_{h}(x)=\tfrac{1}{h}K(\tfrac{x}{h})$, where $K(\cdot)$ is a Kernel function such that $K(x/h)=0$ is $|x|>h$ and $h=h_{\!{}_{T}}$ is the smoothing parameter that tends to zero as $T$ diverges, but slower than $1/T$:

$$h=h_{\!{}_{T}}\to 0\quad\mbox{and}\quad Th_{\!{}_{T}}\to\infty\quad\mbox{as}\,T\to\infty.$$

If we right-multiply (3) by $\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)$ and sum over $t$ we obtain

$\displaystyle\sum_{t=1}^{T}\bm{X}_{t}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)$

$\displaystyle=$

$\displaystyle\sum_{t=1}^{T}\mathbf{A}(\tfrac{t}{T})\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)+\sum_{t=1}^{T}\bm{\varepsilon}_{t}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u).$

Since $\mathbb{E}[\bm{\varepsilon}_{t}\bm{X}_{t-1}^{\top}]=0$,

$$\sum_{t=1}^{T}\mathbf{A}(\tfrac{t}{T})\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)-\mathbf{A}(u)\sum_{t=1}^{T}\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)=\sum_{t=1}^{T}[\mathbf{A}(\tfrac{t}{T})-\mathbf{A}(u)]\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u),$$

	$\displaystyle\|\tfrac{1}{T}\sum_{t=1}^{T}\bm{X}_{t}\bm{X}_{t}^{\top}K_{h}(\tfrac{t}{T}-u)-\mathbf{\Gamma}(0,u)\|=\mathcal{O}_{p}(\tfrac{1}{\sqrt{T\,h}}),$
	$\displaystyle\|\tfrac{1}{T}\sum_{t=1}^{T}\bm{X}_{t}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)-\mathbf{\Gamma}(1,u)\|=\mathcal{O}_{p}(\tfrac{1}{\sqrt{T\,h}}),$

and

	$\displaystyle\|\tfrac{1}{T}\sum_{t=1}^{T}[\mathbf{A}(\tfrac{t}{T})-\mathbf{A}(u)]\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)\|$	$\displaystyle\leq$	$\displaystyle\sup_{t,T}\|\mathbf{A}(\tfrac{t}{T})\|\times\|\tfrac{1}{T}\sum_{t=1}^{T}\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)\|$
		$\displaystyle\leq$	$\displaystyle|\tfrac{t}{T}-u|\,\times\|\mathbf{A}^{(1)}(u)\|\times\|\tfrac{1}{T}\sum_{t=1}^{T}\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)\|$
		$\displaystyle\leq$	$\displaystyle\mathcal{O}(\tfrac{1}{T})\times\mathcal{O}(1)\times\|\mathbf{\Gamma}(0,u)+\mathcal{O}_{p}(\tfrac{1}{\sqrt{T\,h}})\|=\mathcal{O}_{p}(\tfrac{1}{T}),$

it makes sense estimating $\mathbf{A}(u)$ as

$$\widehat{\mathbf{A}}(u)=\Big{[}\sum_{t=1}^{T}\bm{X}_{t}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)\Big{]}\Big{[}\sum_{t=1}^{T}\bm{X}_{t-1}\bm{X}_{t-1}^{\top}K_{h}(\tfrac{t}{T}-u)\Big{]}^{-1}.$$

(4)

If we define

	$\displaystyle\underset{T\times r}{\mathbf{X}_{0}}$	$\displaystyle=$	$\displaystyle\{\underset{r\times 1}{\bm{X}_{0}},\underset{r\times 1}{\bm{X}_{1}},\dots,\underset{r\times 1}{\bm{X}_{T-1}}\}^{\top},$		(5)
	$\displaystyle\underset{T\times r}{\mathbf{X}_{1}}$	$\displaystyle=$	$\displaystyle\{\underset{r\times 1}{\bm{X}_{1}},\underset{r\times 1}{\bm{X}_{2}},\dots,\underset{r\times 1}{\bm{X}_{T}}\}^{\top},\,\,\mbox{and}$		(6)
	$\displaystyle\underset{T\times T}{\mathbf{K}_{\!{}_{T}}(u)}$	$\displaystyle=$	$\displaystyle{\rm{diag}}\{K_{h}(\tfrac{1}{T}-u),K_{h}(\tfrac{2}{T}-u),\dots,K_{h}(1-u)\},$		(7)

the estimator in (4) can be written as

$$\widehat{\mathbf{A}}(u)=[\mathbf{X}_{1}^{\top}\mathbf{K}_{\!{}_{T}}(u)\mathbf{X}_{0}]\,[\mathbf{X}_{0}^{\top}\mathbf{K}_{\!{}_{T}}(u)\mathbf{X}_{0}]^{-1}.$$

(8)

The estimator in (8) can be obtained as the minimizer of a WLS problem. Consider the LS-VAR(1) in (3), and let us now allow for a time-varying (non-zero) mean

$$\bm{X}_{t}-\bm{\mu}(\tfrac{t}{T})=\mathbf{A}(\tfrac{t}{T})[\bm{X}_{t-1}-\bm{\mu}(\tfrac{t-1}{T})]+\bm{\varepsilon}_{t},\quad t=1,\dots,T,$$

(9)

with $\bm{X}_{0}=\bm{\mu}(0)=\bm{0}$ and $\mathbb{E}[\bm{\varepsilon}_{t}\bm{X}_{s}^{\top}]=\mathbbm{1}_{\{s=t\}}\mathbf{\Gamma}_{\varepsilon}$. If the largest eigenvalue of $\mathbf{A}$ lies inside the unit circle

$$\sup_{u\in(0,1)}\left|{\rm v}_{1}[\mathbf{A}(u)]\,\right|<1,$$

(10)

$\bm{X}_{t}$ is locally stationary and causal. Our goal is to estimate $\mathbf{A}(u)$ at a fixed $u\in(0,1)$ by WLS. We can rewrite (9) as

$$\bm{X}_{t}=\bm{m}(\tfrac{t}{T})+\mathbf{A}(\tfrac{t}{T})\bm{X}_{t-1}+\bm{\varepsilon}_{t},$$

(11)

where $\bm{m}(\tfrac{t}{T})=\bm{\mu}(\tfrac{t}{T})-\mathbf{A}(\tfrac{t}{T})\bm{\mu}(\tfrac{t-1}{T})$ for all $t=1,\dots,T$. If we define

	$\displaystyle\underset{r\times(r+1)}{\mathbf{B}(u)}$	$\displaystyle=$	$\displaystyle[\bm{m}(u),\mathbf{A}(u)]\quad\mbox{and}$
	$\displaystyle\underset{(r+1)\times 1}{\bm{Z}_{t}}$	$\displaystyle=$	$\displaystyle\begin{pmatrix}1\\ \bm{X}_{t-1}\end{pmatrix}$

for all $t=1,\dots,T$, model (11) can be written as

$$\bm{X}_{t}=\mathbf{B}(\tfrac{t}{T})\bm{Z}_{t}+\bm{\varepsilon}_{t},$$

(12)

and it makes sense to define $\widehat{\mathbf{B}}(u)$ as the minimizer of the weighted loss function

$$\sum_{t=1}^{T}\|\bm{X}_{t}-\mathbf{B}(\tfrac{t}{T})\bm{Z}_{t}\|^{2}\,K_{h}(\tfrac{t}{T}-u),$$

(13)

where the bandwidth sequence $h\equiv h_{T}$ tends to zero slower than $T^{-1}$: $h_{T}\to 0$ and $T\,h_{T}\to\infty$ as $T\to\infty$. The following proposition provides a closed-form of the local-constant minimizer of (13). We consider model (12) and use the local-constant approximation of $\mathbf{B}(\tfrac{t}{T})$ in a neighborhood of $u$

$$\mathbf{B}(\tfrac{t}{T})\approx\mathbf{B}(u)$$

to estimate $\mathbf{B}(u)$, our parameter of interest, in the approximate model

$$\bm{X}_{t}\approx\mathbf{B}(u)\bm{Z}_{t}+\bm{\varepsilon}_{t}.$$

Our first result generalizes the Yule-Walker solutions in (8) to allow for the time-varying mean vector $\bm{\mu}$.

The following proposition provides a closed-form of the local-linear minimizer of (13). We consider model (12) and use the local-linear approximation of $\mathbf{B}(\tfrac{t}{T})$ in a neighborhood of $u$

$$\mathbf{B}(\tfrac{t}{T})\approx\mathbf{B}(u)+(\tfrac{t}{T}-u)\mathbf{B}^{(1)}(u)$$

(21)

to estimate $\mathbf{B}(u)$, our parameter of interest, in the approximate model

$$\bm{X}_{t}\approx\mathbf{B}(u)\bm{Z}_{t}+(\tfrac{t}{T}-u)\mathbf{B}^{(1)}(u)\bm{Z}_{t}+\bm{\varepsilon}_{t}.$$

(22)

We first consider the zero-mean locally stationary VAR(1) model in (3), and estimate the VAR matrix by means of the localized Yule-Walker equations. Then we consider the locally stationary VAR(1) model in (9) with time-varying mean, and compare the WLS estimates obtained with local-constant and local-linear weights, respectively.

We simulate model (3) with $r=6$. For $j=1,4$ and $k=1,\dots,r$, we generate the time varying entries of the $r\times r$ matrix $\mathbf{A}(\tfrac{t}{T})$ as

	$\displaystyle A_{j,k}(\tfrac{t}{T})$	$\displaystyle=$	$\displaystyle a_{1}\tfrac{\sqrt{j+3}}{\log(k+3)}\sin\big{(}4\pi\tfrac{t}{T}\tfrac{\sqrt{j+4}}{\log(k+4)}\big{)},$
	$\displaystyle A_{j+1,k}(\tfrac{t}{T})$	$\displaystyle=$	$\displaystyle a_{1}\tfrac{\sqrt{j+2}}{\log(k+3)}\cos\big{(}2\pi\tfrac{t}{T}\tfrac{\sqrt{j+4}}{\log(k+2)}\big{)},$		(32)
	$\displaystyle A_{j+2,k}(\tfrac{t}{T})$	$\displaystyle=$	$\displaystyle a_{2}\tfrac{\sqrt{j+1}}{\log(k+3)}\sin\big{(}\pi\tfrac{t}{T}\tfrac{\sqrt{j+4}}{\log(k+2)}\big{)},$

with $T=800$, $a_{1}=0.2$ and $a_{2}=0.1$. This specification is such that, for all $u\in(0,1)$,

$$0.1<\left|{\rm v}_{1}[\mathbf{A}(u)]\,\right|<0.9.$$

We estimate the parameters according to (8), with $h=0.03$ and the Gaussian Kernel $K(x)=\tfrac{1}{\sqrt{2\pi}}\exp(-0.5\,x^{2})$. The results are reported in Figure 1.

We simulate model (9) with $r=3$. For $k=1,\dots,r$, we generate the time varying entries of the $r\times 1$ vector $\bm{\mu}(\tfrac{t}{T})$ as

$$\mu_{k}(\tfrac{t}{T})=\sqrt{6}\,\sin(\pi\omega_{k}\,\tfrac{t}{T}-\phi_{k})$$

(33)

and the $r\times r$ matrix $\mathbf{A}(\tfrac{t}{T})$ as

	$\displaystyle A_{1,k}(\tfrac{t}{T})$	$\displaystyle=$	$\displaystyle a_{1}\tfrac{\sqrt{6}}{\log(k+3)}\sin\big{(}1.2+2\pi\tfrac{t}{T}\tfrac{\sqrt{7}}{\log(k+4)}\big{)},$
	$\displaystyle A_{2,k}(\tfrac{t}{T})$	$\displaystyle=$	$\displaystyle a_{1}\tfrac{\sqrt{5}}{\log(k+3)}\cos\big{(}1.2+2\pi\tfrac{t}{T}\tfrac{\sqrt{7}}{\log(k+2)}\big{)},$		(34)
	$\displaystyle A_{3,k}(\tfrac{t}{T})$	$\displaystyle=$	$\displaystyle a_{2}\tfrac{\sqrt{4}}{\log(k+3)}\sin\big{(}1.2+\pi\tfrac{t}{T}\tfrac{\sqrt{7}}{\log(k+2)}\big{)},$

with $T=600$, $a_{1}=0.3$ and $a_{2}=0.2$, $\omega_{k}=0.5+k$, and $\phi_{k}=0.2+k/3$. This specification is such that, for all $u\in(0,1)$,

$$0<\left|{\rm v}_{1}[\mathbf{A}(u)]\,\right|<0.9.$$

Figure 2 exhibits the parameters estimated by WLS with $h=0.04$ and the Epanechnikov Kernel $K(x)=\tfrac{3}{4}(1-u^{2})\mathbbm{1}_{\{|u|\leq 1\}}(x)$. The local constant estimates, obtained according to (15)-(16), are presented reported in Figure 2(a). The local linear estimates, obtained according to (25)-(27), are reported in Figure 2(b).

Figure 1 shows that in the absence of the (time-varying) mean, that is when $\bm{\mu}(u)\equiv 0$, the local constant estimator performs very well. However, as Figure 2(a) illustrates, this is not the case in the presence of a (time-varying) mean.

It is clear from Figure 2 that although the local constant estimates look satisfactory, the local linear approach delivers superior results. As in the univariate case, the bias of the local linear estimator only depends on the second derivative (of the unknown regression function), and this does not come as a cost in the asymptotic variance. Moreover, the local linear estimator does not suffer from boundary bias problems. The quality of the estimates in Figure 2(b) is remarkable.

In this paper we consider the problem of estimating the time-varying mean-vector and the time-varying AR-matrix of a locally stationary VAR(1). We provide the closed form definition of the local linear solution to the weighted least squares, in a way that the mean-vector and the AR-matrix are estimated jointly.

The asymptotic properties of our estimator need to be studied. Also, it would be interesting to develop data driven methods to select the smoothing parameters. Moreover, we might consider the problem of estimating the parameters of a locally stationary VAR($p$) of order $p>1$.

An important contribution from future studies is the extension of the WLS in (40) to the high-dimensional setting: $r>>T$. To this end, the WLS approach can be can be generalized in more than one direction. In fact, the closed form in (25)-(27) of Theorem 1 becomes particularly attractive in the case the length $r$ of the time series becomes large. Indeed, we can stick to the linear regression model in (22) with the same assumptions as in Section 3, and fit (22) in a way to shrink the regression coefficients towards zero. More precisely we can consider minimizing, with respect to $\mathbf{B}$, the following WLS-Ridge loss-function

$$\|\mathbf{X}_{1}-\widetilde{\mathbf{Z}}_{0}\mathbf{B}(u)^{\top}\|^{2}_{\mathbf{K}(u)}+\lambda\|{\mathbf{B}}(u)\|^{2},$$

(35)

the effect of the penalty being to shrink the entries of $\mathbf{B}$ towards zero. The approach based on (35) can be generalized to the case of a non-spherical penalty. The loss function corresponding to this scenario is

$$\|\mathbf{X}_{1}-\widetilde{\mathbf{Z}}_{0}{\mathbf{B}}(u)^{\top}\|^{2}_{\mathbf{K}(u)}+\|{\mathbf{B}}(u)\|^{2}_{\mathbf{\Lambda}},$$

(36)

which comprises a WLS criterion – as in (40) – and a generalized ridge penalty given by the matrix $\mathbf{\Lambda}(u)$. In both (35) and (36), the $(T\times T)$ matrix $\mathbf{K}(u)$ is diagonal, the $t$-th element $\{K_{t}(u)=\tfrac{1}{h}K(\tfrac{u-t/T}{h}),\,1\leq t\leq T\}$ representing the weight of the $t$-th observation such that $\tfrac{1}{T}K_{t}(u)\in[0,1]$. The penalty in (36) is a quadratic form with penalty parameter $\mathbf{\Lambda}$, an $r$-dimensional positive-definite matrix. When $\mathbf{\Lambda}=\lambda\mathbf{I}_{r}$, we obtain the spherical penalty of the WLS-Ridge regression in (35). Generalizing the (positive) scalar $\lambda$ to the class of (positive definite) matrices $\mathbf{\Lambda}$ allows for (i) different penalization per regression parameter, and (ii) joint shrinkage among the elements of $\mathbf{B}(u)$.