Semi-parametric estimation and prediction intervals in state space models

TL;DR

This paper develops semi-parametric methods for estimating unknown components of state space models and constructing reliable prediction intervals, addressing limitations of traditional parametric approaches.

Contribution

It introduces consistent estimation techniques for unknown state transition matrices and noise distributions, along with an algorithm for prediction intervals in semi-parametric state space models.

Findings

Estimation methods are consistent under unknown parameters.

Constructed prediction intervals are asymptotically valid.

Numerical experiments demonstrate effectiveness of the proposed methods.

Abstract

Literatures in state space models focus on parametric inference and prediction, which fail if the state space model is not fully specified and the maximum likelihood estimation does not work. In this paper, we assume the state transition matrix and the distribution of state noises are unknown. Under this assumption, we provide methods to consistently estimate these terms. In addition, we introduce an algorithm to construct consistent prediction intervals for state vectors and future observations. We complement the asymptotic results with several numerical experiments.

Tables

Table 1. Table 1: Performance of presented estimators, we perform 500 experiments and record the sample mean (and 90 % percent 90 90\% quantile for crucial estimators) of each item. We do not record errors in estimating f X subscript 𝑓 𝑋 f_{X} for (O1) and (O2) since we do not know the distribution of X 1 subscript 𝑋 1 X_{1} when ϵ 1 subscript italic-ϵ 1 \epsilon_{1} is not normal

Model	size	${\Vert \hat{A}-A\Vert }_2$	${\Vert {\hat{f}}_X-f_X\Vert }_{T_2}$	${\Vert {\hat{f}}_{ϵ}-f_{ϵ}\Vert }_{T_2}$		${\Vert \hat{z}-z\Vert }_{T_2}$		${\Vert \hat{g}-g\Vert }_{T_2}$
		mean	mean	mean	$90\%$ quan	mean	$90\%$ quan	mean	$90\%$ quan
(O1)	500	0.040	/	0.153	0.230	0.072	0.112	0.039	0.059
(O1)	2000	0.021	/	0.124	0.174	0.053	0.077	0.030	0.046
(S1)	500	0.043	0.059	0.164	0.223	0.045	0.076	0.021	0.037
(S1)	2000	0.020	0.046	0.142	0.171	0.025	0.042	0.013	0.021
(O2)	5000	0.062	/	0.173	0.211	0.080	0.102	0.050	0.078
(S2)	5000	0.064	0.019	0.041	0.052	0.025	0.034	0.010	0.017

Table 2. Table 2: Performance on prediction, the result is calculated through 500 simulations, nominal coverage probability is 1 − α = 0.95 1 𝛼 0.95 1-\alpha=0.95 , F, PX and PY respectively represents filtering, predicting state vectors and predicting future observations.

Model	size	Method	coverage probability			mean length
			F	PX	PY	F	PX	PY
(O1)	500	Algorithm 1	0.944	0.922	0.904	4.239	4.942	6.113
(O1)	500	Example 2	0.838	0.906	0.920	4.588	5.370	6.649
(O1)	500	Bootstrap	/	/	0.964	/	/	7.550
(O1)	2000	Algorithm 1	0.940	0.934	0.948	4.224	5.403	6.605
(O1)	2000	Example 2	0.880	0.904	0.912	4.588	5.370	6.649
(O1)	2000	Bootstrap	/	/	0.972	/	/	7.602
(S1)	500	Algorithm 1	0.960	0.926	0.910	3.846	4.776	5.738
(S1)	500	Example 2	0.838	0.896	0.906	4.588	5.370	6.649
(S1)	500	Bootstrap	/	/	0.968	/	/	7.433
(S1)	2000	Algorithm 1	0.960	0.924	0.930	3.847	5.050	6.091
(S1)	2000	Example 2	0.842	0.908	0.922	4.588	5.370	6.649
(S1)	2000	Bootstrap	/	/	0.980	/	/	7.519

Equations

$$\text{(S1)}\{\begin{array}{l}{X}_{n+1}={ϵ}_{n+1},{ϵ}_1\sim \mathcal{N}(0,1)\\ Y_n=X_n+{\eta }_n,{\eta }_1\sim \mathcal{N}(0,2)\end{array}\text{(S2)}\{\begin{array}{l}{X}_{n+1}={ϵ}_{n+1},{ϵ}_1\sim \mathcal{N}(0,2)\\ Y_n=X_n+{\eta }_n,{\eta }_1\sim \mathcal{N}(0,1)\end{array}$$

$$\text{(S3)}\{\begin{array}{l}{X}_{n+1}={ϵ}_{n+1},{ϵ}_1\sim \mathcal{N}(0,4)\\ Y_n=\frac{1}{2}{X}_n+{\eta }_n,{\eta }_1\sim \mathcal{N}(0,2)\end{array}\text{(S4)}\{\begin{array}{l}{X}_{n+1}={ϵ}_{n+1},{ϵ}_1\sim \mathcal{N}(0,I_2)\\ Y_n=[1,0]X_{n+1}+{\eta }_n,{\eta }_1\sim \mathcal{N}(0,2)\end{array}$$

$$X_{n+1}=A{X}_n+{ϵ}_{n+1}$$

$$Y_n=B{X}_n+{\eta }_n,n=...,-1,0,1,...$$

$$\mathcal{F}f(x)={\int }_{{\mathbf{R}}^d}f(y)\exp (i{y}^{T}x)dy$$

$${\mathcal{F}}^{-1}f(x)=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbf{R}}^d}f(y)\exp (-i{y}^{T}x)dy$$

$$\sum _{j=1}^d{\int }_{{\mathbf{R}}^d}{x}_j^4{f}_{\tau }(x_1,...,x_d)d{x}_1...d{x}_d<\infty$$

$${\int }_{{\mathbf{R}}^d}|\mathcal{F}{f}_{\tau }(x_1,...x_d)|d{x}_1...d{x}_d<\infty$$

$$K_h(x_1,...,x_d)=\prod _{i=1}^d\frac{1}{h}G\left(\frac{x_i}{h}\right)\Rightarrow \mathcal{F}{K}_h(x_1,...,x_d)=\prod _{i=1}^d\mathcal{F}G(h{x}_i)$$

$$|\mathcal{F}{f}_{\eta }(x_1,...x_d)|\ge c\prod _{j=1}^d(x_j^2+1{)}^{-\beta /2},\sum _{j=1}^d{\int }_{{\mathbf{R}}^d}|\mathcal{F}{f}_{ϵ}(x_1,...,x_d){|}^2(1+x_j^2{)}^{b/2}d{x}_1...d{x}_d<\infty$$

$$|\mathcal{F}{f}_{\eta }(x_1,...,x_d)|\ge c\prod _{j=1}^d\exp (-\gamma |x_j{|}^{\beta }),\sum _{j=1}^d{\int }_{{\mathbf{R}}^d}|\mathcal{F}{f}_{ϵ}(x_1,...,x_d){|}^2\exp (r|x_j{|}^b)d{x}_1...d{x}_d<\infty$$

$$G(x)=\frac{1}{\pi x^2}\left(\cos (x)-\cos (2x)\right)\Rightarrow \mathcal{F}G(x)=\{\begin{array}{l}x+2x\in [-2,-1]\\ 1x\in [-1,1]\\ 2-xx\in [1,2]\\ 0\text{Otherwise}\end{array}$$

$$\mathcal{F}{f}_X(x)=\mathcal{F}{f}_X(A^{T}x)\mathcal{F}{f}_{ϵ}(x)\Rightarrow |\mathcal{F}{f}_X(x)|\le |\mathcal{F}{f}_{ϵ}(x)|$$

$$\mathcal{F}{f}_Y(x)=\mathcal{F}{f}_X(B^{T}x)\mathcal{F}{f}_{\eta }(x),Y_{n+1}-BA{B}^{-1}{Y}_n={\eta }_{n+1}+B{ϵ}_{n+1}-BA{B}^{-1}{\eta }_n$$

$$\Rightarrow \mathcal{F}{f}_{Y_2-BA{B}^{-1}{Y}_1}(x)=\mathcal{F}{f}_{\eta }(x)\mathcal{F}{f}_{ϵ}(B^{T}x)\mathcal{F}{f}_{\eta }(-B^{-1T}{A}^T{B}^{T}x)$$

$$\hat{A}=(\sum _{k=3}^n{B}^{-1}{Y}_k{Y}_{k-2}^T{B}^{-1T})(\sum _{k=3}^n{B}^{-1}{Y}_{k-1}{Y}_{k-2}^T{B}^{-1T}{)}^{+}$$

$${\tilde{f}}_X(x)=\frac{1}{n\times (2\pi {)}^d}\sum _{j=1}^n{\int }_{{\mathbf{R}}^d}\frac{\mathcal{F}{K}_h(y)}{\mathcal{F}{f}_{\eta }(B^{-1T}y)}\exp (i{y}^T{B}^{-1}{Y}_j-i{y}^{T}x)dy,{\hat{f}}_X(x)=\max (\Re ({\tilde{f}}_X(x)),0)$$

$${\tilde{f}}_{ϵ}(x)=\frac{1}{(n-1)\times (2\pi {)}^d}\sum _{j=1}^{n-1}{\int }_{{\mathbf{R}}^d}\frac{\mathcal{F}{K}_h(y)}{\mathcal{F}{f}_{\eta }(B^{-1T}y)\mathcal{F}{f}_{\eta }(-B^{-1T}{\hat{A}}^{T}y)}\exp (i{y}^T{B}^{-1}{Y}_{j+1}-i{y}^T\hat{A}{B}^{-1}{Y}_j-i{y}^{T}x)dy$$

$${\hat{f}}_{ϵ}(x)=\max (\Re ({\tilde{f}}_{ϵ}(x)),0)$$

$$I({\tilde{f}}_X(x))={\left(\frac{a}{h\pi }\right)}^d\frac{1}{nR}\sum _{k=1}^R\sum _{j=1}^n\frac{\mathcal{F}{K}_h({\zeta }_k)}{\mathcal{F}{f}_{\eta }(B^{-1T}{\zeta }_k)}\exp (i{\zeta }_k^T{B}^{-1}{Y}_j-i{\zeta }_k^{T}x)$$

$$I({\tilde{f}}_{ϵ}(x))={\left(\frac{a}{h\pi }\right)}^d\frac{1}{(n-1)R}\sum _{k=1}^R\sum _{j=1}^{n-1}\frac{\mathcal{F}{K}_h({\zeta }_k)\exp (i{\zeta }_k^T{B}^{-1}{Y}_{j+1}-i{\zeta }_k^T\hat{A}{B}^{-1}{Y}_j-i{\zeta }_k^{T}x)}{\mathcal{F}{f}_{\eta }(B^{-1T}{\zeta }_k)\mathcal{F}{f}_{\eta }(-B^{-1T}{\hat{A}}^T{\zeta }_k)}$$

$$\underset{x\in {\mathbf{R}}^d}{\sup }|{\hat{f}}_{\tau }(x)-f_{\tau }(x)|\le \underset{x\in {\mathbf{R}}^d}{\sup }{\int }_{{\mathbf{R}}^d}|{\hat{f}}_{\tau }(y)-f_{\tau }(y)|\times |g(x-y)|dy\le \Vert g{\Vert }_{L_2}\times \Vert {\hat{f}}_{\tau }-f_{\tau }{\Vert }_{L_2}=o_p(1)$$

$$Prob\left(Z_n\in \Gamma \right)\to 1-\alpha \text{as }n\to \infty$$

$$(Y_{n+1}-BA{\hat{X}}_n{)}^T(\Pi +B\Sigma B^T+BA{\Omega }_n{A}^T{B}^T{)}^{-1}(Y_{n+1}-BA{\hat{X}}_n)\le c_{{\chi }_d^2,1-\alpha }$$

$$Prob\left(\Vert X_n-B^{-1}{Y}_n{\Vert }_{\infty }\le x\right)={\int }_{\Vert B^{-1}y{\Vert }_{\infty }\le x}{f}_{\eta }(y)dy$$

$$Prob\left(\Vert X_{n+1}-A{B}^{-1}{Y}_n{\Vert }_{\infty }\le x\right)={\int }_{{\mathbf{R}}^d}\left({\int }_{\Vert z{\Vert }_{\infty }\le x}{f}_{ϵ}(z+A{B}^{-1}y)dz\right)f_{\eta }(y)dy$$

$$Prob\left(\Vert Y_{n+1}-BA{B}^{-1}{Y}_n{\Vert }_{\infty }\le x\right)=\frac{1}{|\det (B)|}{\int }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}\left({\int }_{\Vert y{\Vert }_{\infty }\le x}{f}_{ϵ}(B^{-1}y+A{B}^{-1}w-B^{-1}z)dy\right)f_{\eta }(w)f_{\eta }(z)dwdz$$

$$H(x)=\underset{R\to \infty }{\lim }\frac{1}{R}\sum _{j=1}^R{1}_{\Vert B^{-1}{\zeta }_j{\Vert }_{\infty }\le x}=Prob\left(\Vert B^{-1}{\zeta }_1{\Vert }_{\infty }\le x|{\underline{Y}}_n\right)={\int }_{\Vert B^{-1}y{\Vert }_{\infty }\le x}{f}_{\eta }(y)dy$$

$$N(x)=(2x{)}^d\underset{R\to \infty }{\lim }\frac{1}{R}\sum _{j=1}^R{\hat{f}}_{ϵ}({\kappa }_j+\hat{A}{B}^{-1}{\zeta }_j)=(2x{)}^d\mathbf{E}{\hat{f}}_{ϵ}({\kappa }_1+\hat{A}{B}^{-1}{\zeta }_1)|{\underline{Y}}_n$$

$$={\int }_{\Vert z{\Vert }_{\infty }\le x,y\in {\mathbf{R}}^d}{\hat{f}}_{ϵ}(z+\hat{A}{B}^{-1}y)f_{\eta }(y)dzdy$$

$$G(x)=\frac{(2x{)}^d}{|\det (B)|}\underset{R\to \infty }{\lim }\sum _{j=1}^R{\hat{f}}_{ϵ}(B^{-1}{\kappa }_j+\hat{A}{B}^{-1}{\zeta }_{1,j}-B^{-1}{\zeta }_{2,j})=\frac{(2x{)}^d}{|\det (B)|}\mathbf{E}{\hat{f}}_{ϵ}(B^{-1}{\kappa }_1+\hat{A}{B}^{-1}{\zeta }_{1,1}-B^{-1}{\zeta }_{2,1})|{\underline{Y}}_n$$

$$=\frac{1}{|\det (B)|}{\int }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}\left({\int }_{\Vert y{\Vert }_{\infty }\le x}{\hat{f}}_{ϵ}(B^{-1}y+\hat{A}{B}^{-1}w-B^{-1}z)dy\right)f_{\eta }(w)f_{\eta }(z)dwdz$$

$$\underset{0<x\le x_0}{\sup }|N(x)-Prob\left(\Vert X_2-A{B}^{-1}{Y}_1{\Vert }_{\infty }\le x\right)|=o_p(1)$$

$$\underset{0<x\le x_0}{\sup }|G(x)-Prob\left(\Vert Y_2-BA{B}^{-1}{Y}_1{\Vert }_{\infty }\le x\right)|=o_p(1)$$

$$Prob\left(\Vert X_2-A{B}^{-1}{Y}_1{\Vert }_{\infty }\le x_0\right)>1-\alpha /4\text{and}Prob\left(\Vert Y_2-BA{B}^{-1}{Y}_1{\Vert }_{\infty }\le x_0\right)>1-\alpha /4$$

$$Prob\left(\Vert X_n-B^{-1}{Y}_n{\Vert }_{\infty }\le x_{\mathcal{H}}^{\ast }\right)-(1-\alpha )=0$$

Full text

1 Introduction

State space models (along with the Kalman filter [1]) were originally developed to control linear systems (Davis and Vinter [2]) but have been widely used in many disciplines, including econometrics (Stock and Watson [3] and Hamilton [4]), computer science (Rangapuram et.al. [5]), engineering [6] and etc. For some practical applications, Scott et.al. [7] applied Kalman filter to obtain information from video scenes, Suk.et.al. [8] applied state space model to analyze fMRI and Glickman and Stern [9] developed a model to predict NFL scores.

State space models are powerful tools for statisticians as well. For instance, in time series literature, linear state space model can be used to find the best linear mean-square predictors for the future observations (e.g. chapter 12.2, Brockwell and Davis [10]). There are lots of researches in state space models, including maximum likelihood estimation(chapter 7 in Durbin and Koopman [11] and Jong [12]) and Bayesian estimation (chapter 13 in [13]) of parameters. Classical state space models are linear and have normal noises, but Carter and Kohn [14] introduced how to perform inference on non-Gaussian state space models through Gibbs sampling. Briers et.al. [15] introduced how to smooth a non-Gaussian and non-linear state space model. Advancements in sampling techniques like Pompe et.al. [16] also bring new insights in this field.

Despite great achievements in parametric and Bayesian settings, to the best of our knowledge, few researches on state space models have been conducted under non-parametric settings. However, there are interesting problems in this field. For example, by adopting notations in chapter 1.3.3, [13], if we do not know one or many of the following terms: state transition matrix, measurement transition matrix, distribution of state noise and distribution of measurement noise, we hope to find consistent estimators of them. Moreover, major usages of state space models include estimating state vectors based on observations, so we hope to know prediction intervals for state vectors or future observations apart from predictors in the nonparametric prediction setting introduced in Stein [17], Thombs and Schucany [18], Pan and Politis [19] and Politis [20]. Anderson et.al. [21] provided a consistent estimator for state transition matrix of a special state space model, Costa and Alpuim [22] tried to estimate some parameters in the state space model non-parametrically, and Stoffer and Wall [23] developed a bootstrap algorithm for the maximum-likelihood estimator of state space models. Rodriguez and Ruiz tried to find the prediction interval for future observations through bootstrap [24] but did not provide a theoretical justification on why their algorithm worked.

We are interested in the aforementioned problems, but would like to set constraints on the state space models and explain why we need the constraints at first. Example 1 shows that two different state space models can generate the same observations. Specifically, a state space model can not be identified from observations if 1. distributions of state noise and measurement noise are unknown ((S1) and (S2)), 2. state transition matrix and measurement transition matrix are unknown ((S1) and (S3)) and 3. the dimension of state vectors are larger than the dimension of observations((S1) and (S4)). After excluding these situations, we may find methods to consistently estimate and predict a state space model.

Example 1**.**

Consider the following 4 state space models

[TABLE]

Here $\mathcal{N}(0,a)$ denotes multivariate normal distribution with mean [math] and covariance matrix $a$, and if $a$ is a scalar, then it denotes normal distribution with mean [math] and variance $a$. $I_{2}$ is the two dimensional identity matrix. $\epsilon_{i}$ and $\eta_{i}$, $i=...,-1,0,1,...$ are independent and identically distributed and $\epsilon_{i}$ is independent of $\eta_{j}$ for any $i,j$. These four models are different but their observations $Y_{i},\ i=1,2,...,n$ have the same joint distributions. In other words, a set of observations $Y_{i},\ i=1,2,...,n$ can be generated from (S1) to (S4).

In this paper, we assume that the measurement transition matrix and the distribution of measurement noise are known in addition to a set of observations from a state space model. Based on these information, we will provide a method to consistently estimate the state transition matrix and the distribution of state noise. Moreover, focus on filtering and prediction problems, we will provide consistent prediction intervals for unobserved state vectors and future observations.

The rest of this paper is organized as follows. In chapter 2 we introduce the problem setting and the frequently used notations and assumptions. In chapter 3 we present consistent estimators of the state transition matrix and the distribution of state noise and the related asymptotic results. In chapter 4 we introduce algorithms to construct prediction intervals and explain why these algorithms maintain consistency. We perform finite sample simulations in chapter 5 and make conclusions in chapter 6. Technique details are presented in the appendix.

2 Preliminary

In the following of this paper, we focus on the linear state space model

[TABLE]

According to chapter 1.3.3 in [13], we call $A$ the state transition matrix, $B$ the measurement transition matrix, $\epsilon_{n}$ the state noise and $\eta_{n}$ the measurement noise, $n=...,-1,0,1,...$. $A$ and $B$ are $d\times d$ matrix and $\eta_{n},\epsilon_{n},X_{n}$ and $Y_{n}$ are $d$-dimensional random vectors. We respectively denote $f_{\epsilon},\ f_{\eta}$ and $f_{X}$ as the density of $\epsilon_{1},\ \eta_{1}$ and $X_{1}$. for a function $f:\mathbf{R}^{d}\to\mathbf{C}$, we define its fourier transformation and inverse Fourier transformation as

[TABLE]

In addition, for a vector $x=(x_{1},...,x_{d})^{T}\in\mathbf{R}^{d}$, we denote its Euclidean norm as $\|x\|_{2}=\sqrt{\sum_{i=1}^{d}x^{2}_{i}}$, its infinity norm as $\|x\|_{\infty}=\max_{i=1,2,...,d}|x_{i}|$. In addition, for a given non-singular matrix $M$, we denote the vector norm $\|x\|_{M}=\|Mx\|_{\infty}$, from the equivalence of norms(e.g. Exercise 5.6 in [25]), we can find constants $c_{M},C_{M}>0$ such that $c_{M}\|x\|_{M}\leq\|x\|_{\infty}\leq C_{M}\|x\|_{M}$ for any $x\in\mathbf{R}^{d}$. For a matrix $M$, we define its $2$-norm as $\|M\|_{2}=\sup_{x\in\mathbf{R}^{d},\|x\|_{2}=1}\|Mx\|_{2}$. For a function $f:\mathbf{R}^{d}\to\mathbf{C}$, we define its $L_{2}$ norm as $\|f\|_{L_{2}}=\sqrt{\int_{\mathbf{R}^{d}}|f(x)|^{2}dx}$. Similar as chapter 1.5.1 in [26], for two sequences of real numbers $a_{n}$ and $b_{n}>0$, $n=1,2,...$, we say that $a_{n}=O(b_{n})$ if there exists a constant $C>0$ such that $a_{n}\leq Cb_{n}$ for $\forall n$ and $a_{n}=o(b_{n})$ if $a_{n}/b_{n}\to 0$ as $n\to\infty$. For two sequences of random variables $X_{n}$ and $Y_{n}>0$, $n=1,2,...$, we say that $X_{n}=O_{p}(Y_{n})$ if for any given $\delta>0$, there exists a constant $C_{\delta}>0$ such that $\sup_{n}Prob\left(|X_{n}|\geq C_{\delta}Y_{n}\right)\leq\delta$ and $X_{n}=o_{p}(Y_{n})$ if $X_{n}/Y_{n}\to_{p}0$ as $n\to\infty$. We denote $\to_{p}$ as convergence in probability. Other symbols will be defined before being used. With these notations, we introduce the main assumptions of this paper.

Assumptions

1). For $\tau=\epsilon,\eta$, in model (2) $\tau_{n}\in\mathbf{R}^{d},n=...,-1,0,1,...$ are independent and identically distributed, and the distribution of $\tau_{1}$ is absolutely continuous with respect to Lebesgue measure with density $f_{\tau}$. $\mathbf{E}|\tau_{1,j}|<\infty$ for $j=1,2,...,d$ and $\mathbf{E}\tau=0$. $\epsilon_{n},n=...,-1,0,1,...$ are independent of $\eta_{m},m=...,-1,0,1,...$. $d\times d$ matrix $A$, $B$ and the covariance matrix $\Sigma$ of $\epsilon_{1}$ are non-singular. In addition, we assume that $\|A\|_{2}<1$, the distribution of $X_{n}$ is stationary(like definition 1.3.3 in [10]). We also assume that the measurement transition matrix $B$ and the density of measurement noise $f_{\eta}$ are known in addition to observations $Y_{i},\ i=1,2,...,n$.

2). For $\tau=\epsilon,\eta$, we assume that

[TABLE]

3). $G:\mathbf{R}\to\mathbf{R}$ is a function such that $\int_{\mathbf{R}}|G(x)|dx<\infty$, and $\mathcal{F}G$’s support belongs to $[-a,a]$ for a constant $a>1$, $0\leq\mathcal{F}G\leq 1$ and $\mathcal{F}G=1$ on $[-1,1]$. For each bandwidth $h>0$, we define the kernel function in $\mathbf{R}^{d}$ as

[TABLE]

4). One of the two following conditions happens:

4.a) (Ordinary smooth) There exists constants $b,\beta,c>0$ such that

[TABLE]

and the bandwidth $h=h(n)$ satisfies $h=o(1)$ and $\frac{1}{nh^{4d\beta+d+2}}=o(1)$.

4.b) (Super smooth) There exists constants $c,\beta,\gamma,b,r>0$ such that

[TABLE]

and the bandwidth $h=h(n)$ satisfies $h=o(1)$ and $\frac{1}{h\log^{1/\beta}(n)}=o(1)$.

Classifying distributions through their smoothness is common in the deconvolution literatures such as Fan [27], Fan [28] and Cao [29]. Assumption 4) is similar with conditions used by Comte and Lacour [30] for solving multivariate deconvolution problems. In remark 1, we derive some important corollaries from aforementioned assumptions.

Remark 1**.**

By letting $\tau=\epsilon,\tau$, from theorem 8.22 and 8.26 in [25] and assumption 2), we know that $\mathcal{F}f_{\tau}$ is continuous and vanishes at infinity, and $\mathcal{F}\mathcal{F}^{-1}f_{\tau}=\mathcal{F}^{-1}\mathcal{F}f_{\tau}=f_{\tau}$ almost surely. Thus, without loss of generality, we may assume that $f_{\tau}$ is continuous(which implies uniformly continuous). According to Plancherel’s theorem, $\|\mathcal{F}f_{\epsilon}\|_{L_{2}}=(2\pi)^{d/2}\|f_{\epsilon}\|_{L_{2}}$.

The second thing we would like to emphasize is that $G$ defined in assumption 3) exists. For example, we can set $a=2$ and define function $G$ as

[TABLE]

For any Borel-measurable set $S$ with Lebesgue measure [math], define set $S-Ax=\{y-Ax|y\in S\},\ \forall x\in\mathbf{R}^{d}$, $S-Ax$ has Lebesgue measure [math] as well, so $Prob\left(X_{1}\in S\right)=\mathbf{E}Prob\left(\epsilon_{1}\in S-AX_{0}|X_{0}\right)=0$, which implies that the distribution of $X_{1}$ is absolutely continuous with respect to Lebesgue measure and the density $f_{X}$ is well-defined. From assumption 1), we have for any $x\in\mathbf{R}^{d}$,

[TABLE]

Therefore, assumption 2) and 4) are also valid for $\mathcal{F}f_{X}$.

3 Semi-parametric estimation in state space models

In this chapter we will find estimators for state transition matrix $A$ and the density of state noises $f_{\epsilon}$ as well as the density of state vectors $f_{X}$. From (2), by defining $\mathcal{F}f_{Y}(x)=\mathbf{E}\exp(ix^{T}Y_{1})$,

[TABLE]

Since $\mathcal{F}f_{Y}$ can be consistently estimated from observations, we may apply deconvolution to achieve this goal. We refer Fan [28], Fan [27] and Devroye [31] as an introduction. Masry [32] considered deconvolution problem for dependent data and Comte and Lacour [30] tried to perform deconvolution for random vectors. We present the estimators and summarize the asymptotic results in theorem 1.

Theorem 1**.**

Suppose assumptions 1) to 4),

1. By defining the estimator $\widehat{A}$ as

[TABLE]

$\|\widehat{A}-A\|_{2}=O_{p}(1/\sqrt{n})$, here we denote $A^{+}$ as the pseudo-inverse of matrix $A$.

2. By defining the estimator $\widehat{f}_{X}(x),\ x\in\mathbf{R}^{d}$ as

[TABLE]

with $\Re(x)$ denotes the real part of a complex number $x$, we have $\mathbf{E}\|\widehat{f}_{X}-f_{X}\|_{L_{2}}^{2}=o(1)$.

3. By defining the estimator $\widehat{f}_{\epsilon}(x),\ \forall x\in\mathbf{R}^{d}$ as

[TABLE]

we have $\|\widehat{f}_{\epsilon}-f_{\epsilon}\|_{L_{2}}=o_{p}(1)$

Remark 2**.**

Even though we present estimators $\widetilde{f}_{X}$ and $\widetilde{f}_{\epsilon}$ through multiple integral, in practice we can derive them through Monte Carlo integration. Notice that $\mathcal{F}K_{h}(y)=0$ for $\|y\|_{\infty}>a/h$, in practice for a given positive integer $R$, we can generate $R\times d$ independent and identically distributed random variables $\zeta_{jk},j=1,2,...,R,k=1,2,...,d$ whose marginal distribution is uniform distribution on $[-a/h,a/h]$. Then we define $\zeta_{i}=(\zeta_{i1},...,\zeta_{id})^{T}$ and derive

[TABLE]

from the strong law of large number, when $\widehat{A}$ is non-singular, from the equivalence of norm, we have $\lim_{R\to\infty}I(\widetilde{f}_{X}(x))=\mathbf{E}I(\widetilde{f}_{X}(x))|Y_{1},...,Y_{n}=\widetilde{f}_{X}(x)$ and $\lim_{R\to\infty}I(\widetilde{f}_{\epsilon}(x))=\mathbf{E}I(\widetilde{f}_{\epsilon}(x))|Y_{1},...,Y_{n}=\widetilde{f}_{\epsilon}(x)$ almost surely. Therefore, by selecting sufficiently large $R$, we may use (14), which is easily calculated, to approximate $\widetilde{f}_{X}$ and $\widetilde{f}_{\epsilon}$.

As a special case, if $A=0$ and $B=I_{d}$, the $d$ dimensional identity matrix, then the state space model (2) degenerates to the usual de-convolution model $Y_{n}=\epsilon_{n}+\eta_{n}$, and the estimator $\widetilde{f}_{X}$ coincides with the frequently used estimator in de-convolution setting.

We emphasize an important corollary from theorem 1. Suppose a function $g$ satisfies $\|g\|_{L_{2}}<\infty$, by defining the convolution $f\star g(x)=\int_{\mathbf{R}^{d}}f(y)g(x-y)dy$, for $\tau=X,\epsilon$, from Cauchy’s inequality

[TABLE]

4 Prediction intervals in state space models

Prediction is not a new topic in statistical inference but has seen resurgence with the advent of statistical learning (Hastie et.al.[33]). We refer Geisser [34] for a comprehensive introduction and Politis [20] for a more recent exposition. Apart from a precise predictor, statisticians always use prediction intervals to quantify uncertainty in a prediction [35]. In this chapter, we will provide an algorithm to find the prediction intervals for state vectors and future observations and theoretically justify why the presented algorithm generates consistent prediction intervals. We adopt definition 2.4.1 in [20] to define a consistent prediction interval.

Definition 1** (Consistent prediction interval).**

Suppose there are $n$ observations $Y_{1},...,Y_{n}$ and the set $\Gamma=\Gamma(Y_{1},...,Y_{n})$ as well as the new random variable $Z_{n}$ satisfy

[TABLE]

then we call the set $\Gamma$ the $1-\alpha$ consistent prediction interval for the random variable $Z_{n}$. Here $0<\alpha<1$.

In this paper, we are interested in $Z_{n}=X_{n}$(filtering) and $Z_{n}=X_{n+1},\ Y_{n+1}$(prediction).

Example 2** (Consistent prediction intervals for standard state space model).**

We adopt the definition of a standard state space model in chapter 1.3.3, [13], i.e. $A,\ B$ and the covariance matrix $\Sigma$ of $\epsilon_{1}$ as well as the covariance matrix $\Pi$ of $\eta_{1}$ are known. $\epsilon_{1},\ \eta_{1}$ are multivariate normal distributed. In addition, we assume $X_{1}=[0,0,...,0]^{T}$ and $Y_{0}=[0,0,...,0]^{T}$. From Kalman prediction (proposition 12.2.2 in [10]), by defining $\widehat{X}_{n}$ as the Kalman one-step predictor for $X_{n}$, we may derive the error covariance matrix $\Omega_{n}=\mathbf{E}(X_{n}-\widehat{X}_{n})(X_{n}-\widehat{X}_{n})^{T}$. Since $\mathbf{E}X_{n}=0$ and $\mathbf{E}\widehat{X}_{n}=0$ (we can prove this result by $\mathbf{E}\widehat{X}_{1}=0$ and induction), if $\Omega_{n}$ is non-singular, one of the $1-\alpha$ consistent prediction interval for $X_{n}$ can be given by $(X_{n}-\widehat{X}_{n})^{T}\Omega_{n}^{-1}(X_{n}-\widehat{X}_{n})\leq c_{\chi^{2}_{d},1-\alpha}$, here $c_{\chi^{2}_{d},1-\alpha}$ is the $1-\alpha$ quantile of $\chi^{2}$ distribution with degree of freedom $d$. Since $X_{n+1}-A\widehat{X}_{n}=A(X_{n}-\widehat{X}_{n})+\epsilon_{n+1}$, the $1-\alpha$ consistent prediction interval for $X_{n+1}$ can be given by $(X_{n+1}-A\widehat{X}_{n})^{T}(A\Omega_{n}A^{T}+\Sigma)^{-1}(X_{n+1}-A\widehat{X}_{n})\leq c_{\chi^{2}_{d},1-\alpha}$. Since $Y_{n+1}-BA\widehat{X}_{n}=\eta_{n+1}+B\epsilon_{n+1}+BA(X_{n}-\widehat{X}_{n})$, the $1-\alpha$ consistent prediction interval for $Y_{n+1}$ can be derived by

[TABLE]

Consistency of the aforementioned prediction intervals relies on two facts: 1. the distribution of $X_{n}-\widehat{X}_{n}$ only depends on the covariance matrix $\Omega_{n}$, and 2. linear combinations of random variables preserve normality. These two facts fail if $\epsilon_{1}$ and $\eta_{1}$ are not normally distributed.

Kalman recursion (chapter 12.2 in [10]) is one of the most popular methods to deal with state space models. However, when the state noise and the measurement noise are not normal, distribution of the predictor will be hard to derive. Cappé et.al. [13] introduced how to filter and smooth a general state space model but they assumed nothing needed to be estimated. We use another method to derive predictors. If $\eta_{i}=0$, then model 2 degenerates to a vector autoregressive model. In this situation, according to Kim [36], predictors for $X_{n}$, $X_{n+1}$ and $Y_{n+1}$ can respectively be $B^{-1}Y_{n}$, $\widehat{A}B^{-1}Y_{n}$ and $B\widehat{A}B^{-1}Y_{n}$. When $\eta_{i}\neq 0$, we still can use these predictors but the distributions of predictive roots(definition 2.4.2 in [20]) are different from vector autogressive models.

Lemma 1**.**

Suppose $Y_{i},i=1,2,...,n$ are generated from state space model 2 and assumption 1) is satisfied, then for $\forall 0<x\in\mathbf{R}$

1. Filtering:

[TABLE]

2. Predicting the state vector:

[TABLE]

3. Predicting the future observation:

[TABLE]

Algorithm 1 performs binary searches on the estimated cumulative distribution functions of predictive roots $X_{n}-B^{-1}Y_{n}$, $X_{n+1}-\widehat{A}B^{-1}Y_{n}$ and $Y_{n+1}-B\widehat{A}B^{-1}Y_{n}$ to find $1-\alpha$ quantiles. Validness of algorithm 1 can be assured if the estimated distributions converge to the distributions of predictive roots. We prove this result in theorem 2 and corollary 1.

Remark 3**.**

From the strong law of large number, by defining $\underline{Y}_{n}=(Y_{1},...,Y_{n})$, almost surely for any given $x\in\mathbf{R}^{d}$,

[TABLE]

which are the desired estimators.

Compare to classical de-convolution literatures like [27] or [30], algorithm 1 avoids calculating multidimensional integrations, which is hard to manipulate. There are several methods to generate random vectors with desired density $f_{\eta}$, e.g., Metropolis-Hastings algorithm [37] and Accept-Reject algorithm [13].

From (13), we know that $\widehat{f}_{\epsilon}$ is bounded, so from dominated convergence theorem $H,N$ and $G$ are continuous and increasing. By defining $M=H,N,G$, if $M(0)<1-\alpha$ and $\lim_{x\to\infty}M(x)>1-\alpha$, from intermediate value theorem there exists $x^{*}_{\mathcal{M}}$ such that $M(x^{*}_{\mathcal{M}})=1-\alpha$. Since $x^{*}_{\mathcal{M}}$ can be approximated in algorithm 1 through setting $R\to\infty$ and $\varepsilon\to 0$, we only consider $M$ and $x^{*}_{\mathcal{M}}$ in theoretical justifications.

Theorem 2**.**

Suppose assumptions 1) to 4), then for arbitrary given $x_{0}>0$

[TABLE]

Here $H,\ N$ and $G$ coincide with remark 3.

For a given $0<\alpha<1$, by choosing sufficiently large $x_{0}$ such that

[TABLE]

and apply theorem 2 on $(0,x_{0}]$, the probability of $x_{\mathcal{N}}^{*}$ and $x_{\mathcal{G}}^{*}$ exist converges to $1$ as the sample size increases. In corollary 1 we explain why $x_{\mathcal{M}}^{*},\ \mathcal{M}=\mathcal{H},\mathcal{N},\mathcal{G}$, the asymptotic values of $x_{\mathcal{M}}$ in algorithm 1, contribute to the consistent prediction intervals.

Corollary 1**.**

Suppose assumptions 1) to 4), for arbitrary given $0<\alpha<1$, by defining $x_{\mathcal{M}}^{*},\mathcal{M}=\mathcal{H},\mathcal{N},\mathcal{G}$ as in remark 3

[TABLE]

5 Simulations and Numerical experiments

In this chapter, we justify theoretical results in chapter 3 and 4 through numerical experiments. We separate discussions into ordinary smooth and super smooth situations and design the following four models:

(O1) $X_{n},Y_{n}$ satisfy $X_{n+1}=0.8X_{n}+\epsilon_{n+1}$ and $Y_{n}=1.0X_{n}+\eta_{n}$. $\epsilon_{1}=\epsilon_{1,1}-\epsilon_{1,2}$, $\epsilon_{1,i},\ i=1,2$ have gamma distribution with shape $3/2$ and scale $1/\sqrt{3}$, $\epsilon_{1,1}$ is independent of $\epsilon_{1,2}$. $\eta_{1}=\eta_{1,1}-\eta_{1,2}$, $\eta_{1,i}.\ i=1,2$ have gamma distribution with shape $1/2$ and scale $1.0$. $\eta_{1,1}$ is independent of $\eta_{1,2}$. In this situation, $\eta_{1}$ and $\epsilon_{1}$ have variance $1$ and $\beta=1$. According to assumption 4.a), we choose the bandwidth $h=1.0/n^{1/8}$

(S1) We choose the same $A$ and $B$ as (O1) but $\epsilon_{1}$ and $\eta_{1}$ have normal distribution with mean [math] and variance $1$. According to assumption 4.b), since $\beta=2$, we choose $h=1/log^{0.1}(n)$.

(O2) $X_{n},Y_{n}$ satisfy

[TABLE]

$\epsilon_{i,1}$ is independent of $\epsilon_{i,2}$, $\eta_{i,1}$ is independent of $\eta_{i,2}$, and $\epsilon_{i,j}$, $\eta_{i,j},\ i=...,-1,0,1,...,\ j=1,2$ respectively has the same distribution as $\epsilon_{1}$ and $\eta_{1}$ in (O1). Since $\beta=1$, we choose $h=1.0/n^{1/8}$.

(S2) We choose the same $A$ and $B$ as (O2) but $\epsilon_{i,j},\ j=1,2$, $\eta_{i,j},\ j=1,2$ have normal distribution with mean $0.0$ and variance $1.0$. Since $\beta=2.0$, we choose $h=1.0/\log^{0.1}(n)$.

In the first experiment, we focus on evaluating the performance of presented estimators. Apart from $\widehat{A},\widehat{f}_{X}$ and $\widehat{f}_{\epsilon}$, we also consider the estimator $\widehat{z}(x)=\int_{y\in\mathbf{R}^{d}}\widehat{f}_{\epsilon}(x+\widehat{A}B^{-1}y)f_{\eta}(y)dy$ on $z(x)=\int_{y\in\mathbf{R}^{d}}f_{\epsilon}(x+AB^{-1}y)f_{\eta}(y)dy$ and $\widehat{g}(x)=\int_{\mathbf{R}^{d}\times\mathbf{R}^{d}}\widehat{f}_{\epsilon}(B^{-1}x+\widehat{A}B^{-1}w-B^{-1}z)f_{\eta}(w)f_{\eta}(z)dwdz$ on $g(x)=\int_{\mathbf{R}^{d}\times\mathbf{R}^{d}}f_{\epsilon}(B^{-1}x+AB^{-1}w-B^{-1}z)f_{\eta}(w)f_{\eta}(z)dwdz$. According to lemma 1, if $\widehat{z}$ is close to $z$ and $\widehat{g}$ is close to $g$, algorithm 1 should be able to find the consistent prediction interval for $X_{n+1}$ and $Y_{n+1}$.

The metric we use in the first experiment is $\|f\|_{T_{2}}=\sqrt{\int_{[-1/h,1/h]^{d}}|f(x)|^{2}dx}$, which approximates $\|f\|_{L_{2}}$ and easily derived through Monte Carlo integrations. In table 1 and figure 1, we observe that $\widehat{f}_{\epsilon}$ does not have high convergence rate. Fortunately, after convoluted with $f_{\eta}$, the estimators $\widehat{z}$ and $\widehat{g}$ have better performance.

We demonstrate the performance of algorithm 1 in table 2. Two alternatives are example 2 and the bootstrap algorithm used in [24]. We do not assume knowing the state transition matrix and the distribution of state noise in algorithm 1 but assume everything about the state space model is known in the two alternatives. The coverage probability of prediction intervals constructed by algorithm 1 approximates the nominal coverage probability $1-\alpha$. Algorithm 1 and example 2 tend to find prediction intervals with coverage probability less than $1-\alpha$ while the bootstrap algorithm [24] tends to find a wide prediction interval.

6 Conclusions

Literatures in state space models focus on parametric estimation and filtering. However, if the state space model is not fully specified and maximum likelihood estimation fails, then methods like Kalman filter cannot be applied to make predictions. This paper brings a method to consistently estimate the state transition matrix and the distribution of state noises when they are unknown. In addition, this paper introduces an algorithm to construct prediction intervals of state vectors and future observations with unknown state transition matrix and distribution of state noises. In practice, fully specified state space model is too restrictive to model real life data, so methods introduced in this paper can be good alternatives to parametric estimation and filtering techniques.

Appendix A Proof of theorem 1

We first prove a lemma.

Lemma 2**.**

Suppose $X_{i},Y_{i},\ i=...,-1,0,1,...$ are two stationary $d$ dimensional random vectors and $X_{i},i=...,-1,0,1,...$ are independent of $Y_{i},i=...,-1,0,1,...$. Define $Z_{i}=X_{i}+Y_{i}$ for $\forall i$, suppose the $\alpha-$mixing coefficients $\alpha_{Z}(k),k\geq 0\in\mathbf{Z}$ of $Z_{i}$ satisfy $\sum_{k=0}^{\infty}\alpha_{Z}(k)<\infty$, then we have

[TABLE]

We adopt the definition in chapter 1.1 of [38] and define the $\alpha-$mixing coefficients of $Z_{i}$ as

[TABLE]

Here $\sigma(S)$ denotes the $\sigma$-field generated by $S$.

Proof of lemma 26.

We define $\overline{\mathcal{H}}_{k}$ as the $\sigma$-field generated by $Z_{k},Z_{k+1},...$ and $\underline{\mathcal{H}}_{k}$ as the $\sigma$-field generated by $...,Z_{k-1},Z_{k}$. For any $x\in\mathbf{R}^{d}$, define $\overline{\mathcal{G}}_{k}$ as the $\sigma$-field generated by $\sin(x^{T}Z_{i})\ \text{and}\ \cos(x^{T}Z_{i}),\ i=k,k+1,...$ and $\underline{\mathcal{G}}_{k}$ as the $\sigma$-field generated by $\sin(x^{T}Z_{i})\ \text{and}\ \cos(x^{T}Z_{i}),\ i=...,k-1,k$. Since $\sin,\cos$ are continuous function, we have $\underline{\mathcal{G}}_{k}\subset\underline{\mathcal{H}}_{k}$ and $\overline{\mathcal{G}}_{k}\subset\overline{\mathcal{H}}_{k}$, correspondingly by defining

[TABLE]

We have $\alpha_{\mathcal{G}}(k)\leq\alpha_{Z}(k)$ for any $k=0,1,...$. Notice that for any $k,j$, from lemma 3 in [38] we have

[TABLE]

Therefore, we have

[TABLE]

Since $x\in\mathbf{R}^{d}$ is arbitrarily chosen, we prove (26). ∎

Now we start proving theorem 1.

proof of theorem 1.

We first prove (11). Define $Z_{k}=B^{-1}Y_{k}$ and $\tau_{k}=B^{-1}\eta_{k}$ for $k=1,2,...,n$, we have $Z_{k}=X_{k}+\tau_{k}$. From assumption 1) and 3), we have $X_{k}=\sum_{i=0}^{\infty}A^{i}\epsilon_{k-i}$, and correspondingly

[TABLE]

This implies $\Sigma_{X}=\mathbf{E}X_{1}X_{1}^{T}$ exists and $\Sigma_{X}=\sum_{i=0}^{\infty}A^{i}\Sigma A^{iT}$ is positive definite. Correspondingly

[TABLE]

For a matrix $A$, we define its Frobenius norm as $\|A\|_{F}=\sqrt{Tr(AA^{T})}$, since

[TABLE]

Similarly we have

[TABLE]

Therefore, from (32) to (34), we know that $\|\frac{1}{n-2}\sum_{k=3}^{n}Z_{k}Z_{k-2}^{T}-A^{2}\Sigma_{X}\|_{2}=O_{p}(1/\sqrt{n})$ and similarly $\|\frac{1}{n-2}\sum_{k=3}^{n}Z_{k-1}Z_{k-2}^{T}-A\Sigma_{X}\|_{2}=O_{p}(1/\sqrt{n})$. Apply corollary 6.3.4 and (5.8.4) in [39], since $A\Sigma_{X}$ is non-singular, for sufficiently large $n$, with large probability $\frac{1}{n-2}\sum_{k=3}^{n}Z_{k-1}Z_{k-2}^{T}$ is invertible, and we have $\|(\frac{1}{n-2}\sum_{k=3}^{n}Z_{k-1}Z_{k-2}^{T})^{-1}-(A\Sigma_{X})^{-1}\|_{2}=O_{p}(1/\sqrt{n})$ as well. Therefore, we have

[TABLE]

Which is of $O_{p}(1/\sqrt{n})$ and we prove the first result.

From Plancherel theorem and assumption 3), we have

[TABLE]

If assumption 4.a) happens, we have for any $i=1,2,...,d$,

[TABLE]

If assumption 4.b) happens, we have for any $i=1,2,...,d$,

[TABLE]

Notice that

[TABLE]

is an ARMA model and $A$ is non-singular with $\|A\|_{2}<1$, so the polynomial

[TABLE]

has all roots with absolute value greater than 1. According to theorem 1 in [40], $\left[Z_{i},\tau_{i}\right]^{T}$ is a mixing process and $\exists 0<s<1$ such that the $\alpha-$mixing coefficients $\alpha_{[Z_{i},\tau_{i}]^{T}}(k)=O(s^{|k|})$ for any $k\in\mathbf{Z}$. Combine with lemma 26 and assumption 3), we have

[TABLE]

Since $B$ is non-singular, by change of variable theorem,

[TABLE]

If assumption 4.a) happens, according to lemma 1. in [30] and (37),

[TABLE]

which tends to [math] as $n\to\infty$. If assumption 4.b) happens, similarly we have

[TABLE]

$h=o(1)\Rightarrow\exp\left(-\frac{r}{h^{b}}\right)=o(1)$. On the other hand, let $z=h\log^{1/\beta}(n)\to\infty$, for sufficiently large $n$,

[TABLE]

which tends to $-\infty$ as $n\to\infty$ and we prove 2.

We define

[TABLE]

Similar with (36) to (38) and (41) to (44), if assumption 4.a) happens,

[TABLE]

Since $Z_{j}-AZ_{j-1}$ is a continuous function of $Z_{j},Z_{j-1}$, it is a measurable function in the $\sigma$-field generated by $...,Z_{j-1},Z_{j}$ or $Z_{j-1},Z_{j},...$, this implies for any $k>1$, the $\alpha$-mixing coefficients satisfy $\alpha_{Z_{j}-AZ_{j-1}}(k)\leq\alpha_{Z}(k-1)\Rightarrow\sum_{k=0}^{\infty}\alpha_{Z_{j}-AZ_{j-1}}(k)<\infty$, correspondingly from lemma 26 and Cauchy inequality

[TABLE]

If assumption 4.b) happens, we define $c_{1}=c_{B^{-1T}}$ and $c_{2}=c_{B^{-1T}A^{T}}$,

[TABLE]

Similar with (45), we have $\mathbf{E}\|\widetilde{f}_{\epsilon}(x)^{\dagger}-f_{\epsilon}(x)\|^{2}_{L_{2}}\to 0$ as $n\to\infty$.

Since for any $x,y\in\mathbf{R}^{d}$,

[TABLE]

from Cauchy’s inequality and theorem 8.22 in [25], by defining $C=(\int_{\mathbf{R}^{d}}\|t\|_{2}f_{\eta}(t)dt)^{2}<\infty$,

[TABLE]

Since $\mathbf{E}\frac{1}{n-1}\sum_{j=1}^{n-1}\|Z_{j}\|_{2}^{2}=\mathbf{E}\|Z_{1}\|_{2}^{2}<\infty\Rightarrow\frac{1}{n-1}\sum_{j=1}^{n-1}\|Z_{j}\|_{2}^{2}=O_{p}(1)$, if 4.a) is satisfied

[TABLE]

combine with (48) and (51), we have

[TABLE]

On the other hand, if 4.b) is satisfied, by letting $s=h\log^{1/\beta}(n)\to\infty$, if $\|\widehat{A}\|_{2}\leq 2^{1/\beta}\|A\|_{2}$, for sufficient large $n$, $\frac{4d\gamma\|B^{-1}\|_{2}^{\beta}\|A\|_{2}^{\beta}\times d^{\beta/2}a^{\beta}}{s^{\beta}}<1/4$, correspondingly

[TABLE]

combine with (49), since $\|\widehat{A}-A\|_{2}\to_{p}0$ and (45), we have

[TABLE]

Thus, we have

[TABLE]

∎

Appendix B Proofs for theorems in chapter 4

Proof of lemma 1.

For the first result, notice that

[TABLE]

For the second result, we have

[TABLE]

For the third result, from the change of variable theorem

[TABLE]

∎

Proof of theorem 2.

Notice that for any given $x_{0}>0$

[TABLE]

From remark 1, for any $\delta>0$, there exists $M>0$ such that $f_{\epsilon}(x)<\delta$ for any $\|x\|_{2}>M$, from theorem 1, there exists $C_{\delta}>0$ such that $\|\widehat{A}-A\|_{2}\leq C_{\delta}/\sqrt{n}$ with probability at least $1-\delta$, if this happens, for sufficiently large $n$

[TABLE]

and for $\forall z,y$ such that $\|z\|_{\infty}\leq x_{0}$ and $\|y\|_{2}>n^{1/4}$, we have $\|z+AB^{-1}y\|_{2}\geq\frac{\|y\|_{\infty}}{C_{AB^{-1}}}-\|z\|_{2}\geq\frac{n^{1/4}}{C_{AB^{-1}}\sqrt{d}}-\sqrt{d}x_{0}>M$ and $\|z+\widehat{A}B^{-1}y\|_{2}\geq\|z+AB^{-1}y\|_{2}-\|\widehat{A}-A\|_{2}\|B^{-1}\|_{2}\|y\|_{2}\geq\left(\frac{1}{C_{AB^{-1}}\sqrt{d}}-\frac{C_{\delta}\|B^{-1}\|_{2}}{\sqrt{n}}\right)\|y\|_{2}-\|z\|_{2}>M$, this implies

[TABLE]

Combine with theorem 1 and (60) we prove the first result.

Notice that for any $\delta>0$, there exists $M>0$ such that $\int_{\|w\|_{\infty}\geq M\cup\|z\|_{\infty}\geq M}f_{\eta}(w)f_{\eta}(z)dwdz<\delta$, then

[TABLE]

If $\|\widehat{A}-A\|_{2}\leq C_{\delta}/\sqrt{n}$, then from remark 1, for sufficiently large $n$,

[TABLE]

which implies the second result. ∎

Proof of corollary 1.

Lemma 1 and remark 3 imply the first result.

Notice that $Y_{n}=O_{p}(1)$, from theorem 1 we know that $\|(\widehat{A}-A)B^{-1}Y_{n}\|_{2}\leq\|\widehat{A}-A\|_{2}\|B^{-1}\|_{2}\|Y_{n}\|_{2}=O_{p}(1/\sqrt{n})$ and $\|B(\widehat{A}-A)B^{-1}Y_{n}\|_{2}\leq\|B\|_{2}\|\widehat{A}-A\|_{2}\|B^{-1}\|_{2}\|Y_{n}\|_{2}=O_{p}(1/\sqrt{n})$, so for any $\delta>0$, there exists a constant $C_{\delta}$ such that $Prob\left(\|(\widehat{A}-A)B^{-1}Y_{n}\|_{2}\leq\frac{C_{\delta}}{\sqrt{n}}\right)>1-\delta$ and $Prob\left(\|B(\widehat{A}-A)B^{-1}Y_{n}\|_{2}\leq\frac{C_{\delta}}{\sqrt{n}}\right)>1-\delta$ for any $n$. We adopt (1.2) in [41] for and define

[TABLE]

for all $0<\alpha<1$. From lemma 1 and dominated convergence theorem we know that $Prob\left(\|X_{2}-AB^{-1}Y_{1}\|_{\infty}\leq x\right)$ and $Prob\left(\|Y_{2}-BAB^{-1}Y_{1}\|_{\infty}\leq x\right)$ are continuous in $x>0$, so $Prob\left(\|X_{2}-AB^{-1}Y_{1}\|_{\infty}\leq c_{\mathcal{N},1-\alpha}\right)=1-\alpha$ and $Prob\left(\|Y_{2}-BAB^{-1}Y_{1}\|_{\infty}\leq c_{\mathcal{G},1-\alpha}\right)=1-\alpha$.

For a given $0<\alpha<1$, by choosing sufficiently large $x_{0}$ such that (23) happens and adopting notations in remark 3, for any $0<\delta<\alpha/8$ such that $1-\alpha-2\delta>0$, (22) implies

[TABLE]

In (66) we consider $Prob\left(\|X_{n+1}-AB^{-1}Y_{n}\|_{\infty}\leq x\right)$ as a function in $x$ and plug in $x=x^{*}_{\mathcal{N}}$ instead of taking expectations. Similarly we can derive $x^{*}_{\mathcal{G}}\geq c_{\mathcal{G},1-\alpha-\delta}$ and $x^{*}_{\mathcal{G}}\leq c_{\mathcal{G},1-\alpha+\delta}$. From theorem 1, we have for sufficiently large $n$

[TABLE]

which implies (24) by choosing $\delta\to 0$. Similarly we can replace $X_{n+1}-\widehat{A}B^{-1}Y_{n}$ and $x_{\mathcal{N}}^{*}$ in (67) by $Y_{n+1}-B\widehat{A}B^{-1}Y_{n}$ and $x^{*}_{\mathcal{G}}$ and derive the third result in (24). ∎