Maximum Likelihood Estimation of the VAR(1) Model Parameters with Missing Observations

Abstract

Missing-data problems are extremely common in practice. To achieve reliable inferential results, we need to take into account this feature of the data. Suppose that the univariate data set under analysis has missing observations. This paper examines the impact of selecting an auxiliary complete data set—whose underlying stochastic process is to some extent interdependent with the former—to improve the efficiency of the estimators for the relevant parameters of the model. The Vector AutoRegressive (VAR) Model has revealed to be an extremely useful tool in capturing the dynamics of bivariate time series. We propose maximum likelihood estimators for the parameters of the VAR(1) Model based on monotone missing data pattern. Estimators’ precision is also derived. Afterwards, we compare the bivariate modelling scheme with its univariate counterpart. More precisely, the univariate data set with missing observations will be modelled by an AutoRegressive Moving Average (ARMA(2,1)) Model. We will also analyse the behaviour of the AutoRegressive Model of order one, AR(1), due to its practical importance. We focus on the mean value of the main stochastic process. By simulation studies, we conclude that the estimator based on the VAR(1) Model is preferable to those derived from the univariate context.

1. Introduction

Statistical analyses of data sets with missing observations have long been addressed in the literature. For instance, Morrison [1] deduced the maximum likelihood estimators of the parameters of the multinormal mean vector and covariance matrix for the monotonic pattern with only a single incomplete variate. The exact expectations and variances of the estimators were also deduced. Dahiya and Korwar [2] obtained the maximum likelihood estimators for a bivariate normal distribution with missing data. They focused on estimating the correlation coefficient as well as the difference of the two means. Following this line of research and having in mind that the majority of the empirical studies are characterised by temporal dependence between observations, we will try to generalise the previous study by introducing a bivariate time series model to describe the relationship between the processes under consideration.

The literature on missing data has expanded in the last decades focusing mainly on univariate time series models [3–7], but there is still a lack of developments in the vectorial context.

This paper aims at analysing the main properties of the estimators from data generated by one of the most influential models in empirical studies, that is, the first-order Vector AutoRegressive (VAR(1)) Model, when the data set from the main stochastic process, designated by {𝑌_𝑡}_𝑡∈ℤ, has missing observations. Therefore, we assume that there is also available a suitable auxiliary stochastic process, denoted by {𝑋_𝑡}_𝑡∈ℤ, which is to some extent interdependent with the main stochastic process. Additionally, the data set obtained from this process is complete. In this context, a natural question arises: is it possible to exchange information between the two data sets to increase knowledge about the process whose data set has missing observations, or should we analyse the univariate stochastic process by itself? The goal of this paper is to answer this question.

Throughout this paper, we assume that the incomplete data set has a monotone missing data pattern. We follow a likelihood-based approach to estimate the parameters of the model. It is worth pointing out that, in the literature, likelihood-based estimation is largely used to manage the problem of missing data [3, 8, 9]. The precision of the maximum likelihood estimators is also derived.

In order to answer the question raised above, we must verify if the introduction of an auxiliary variable for estimating the parameters of the model increases the accuracy of the estimators. To accomplish this goal, we compare the precision of the estimators just cited with those obtained from modelling the dynamics of the univariate stochastic process {𝑌_𝑡}_𝑡∈ℤ by an AutoRegressive Moving Average (ARMA(2,1)) Model, which corresponds to the marginal model of the bivariate VAR(1) Model [10, 11]. The behaviour of the AutoRegressive Model of order one, AR(1), is also analysed due to its practical importance in time series modelling. Simulation studies allow us to assess the relative efficiency of the different approaches. Special attention is paid to the estimator for the mean value of the stochastic process about which information available is scarce. This is a reasonable choice given the importance of the mean function of a stochastic process in understanding the behaviour of the time series under consideration.

The paper is organised as follows. In Section 2, we review the VAR(1) Model and highlight a few statistical properties that will be used in the remaining sections. In Section 3, we establish the monotone pattern of missing data and factorise the likelihood function of the VAR(1) Model. The maximum likelihood estimators of the parameters are obtained in Section 4. Their precision is also deduced. Section 5 reports the simulation studies in evaluating different approaches to estimate the mean value of the stochastic process {𝑌_𝑡}_𝑡∈ℤ. The main conclusions are summarised in Section 6.

2. Brief Description of the VAR(1) Model

In this section, a few properties of the Vectorial Autoregressive Model of order one are analysed. These features will play an important role in determining the estimators for the parameters when there are missing observations, as we will see in Section 4.

Hereafter, the stochastic process underlying the complete data set is denoted by {𝑋_𝑡}_𝑡∈ℤ, while the other one is represented by {𝑌_𝑡}_𝑡∈ℤ. The VAR(1) Model under consideration takes the form 𝑋_𝑡=𝛼₀+𝛼₁𝑋_𝑡−1+𝜖_𝑡,𝑌_𝑡=𝛽₀+𝛽₁𝑌_𝑡−1+𝛽₂𝑋_𝑡−1+𝜉_𝑡,𝑡=0,±1,±2,…, (1) where 𝜖_𝑡 and 𝜉_𝑡 are Gaussian white noise processes with zero mean and variances 𝜎²_𝜖 and 𝜎²_𝜉, respectively. The structure of correlation between the error terms is different from zero only at the same date 𝑡, that is, Cov(𝜖_𝑡−𝑖,𝜉_𝑡−𝑗)=𝜎_𝜖𝜉, for 𝑖=𝑗; Cov(𝜖_𝑡−𝑖,𝜉_𝑡−𝑗)=0, for 𝑖≠𝑗,𝑖,𝑗∈ℤ. Exchanging information between both time series might introduce some noise in the overall process. Therefore, transfer of information from the smallest series to the largest one is not allowed here.

We have to introduce the restrictions |𝛼₁|<1 and |𝛽₁|<1. They ensure not only that the underlying processes are ergodic for the respective means but also that the stochastic processes are covariance stationary (see Nunes [12, ch.3]). Hereafter, we assume that these restrictions are satisfied.

Next, we overview some relevant properties of the VAR(1) Model (1). Theoretical details can be found in Nunes [12, ch.3].

The mean values of 𝑋_𝑡 and 𝑌_𝑡 are, respectively, given by 𝐸𝑋_𝑡=𝛼₀1−𝛼₁,𝐸𝑌_𝑡=𝛼₀𝛽₂+𝛽₀1−𝛼₁1−𝛼₁1−𝛽₁. (2)

Concerning the covariance structure of the process 𝑋_𝑡, Cov𝑋_𝑡−𝑖,𝑋_𝑡−𝑗=𝜎²_𝜖𝛼₁^{|𝑖−𝑗|}1−𝛼²₁,∀_{𝑖,𝑗∈ℤ}. (3)

For 𝛼₁≠𝛽₁, the covariance of the stochastic process 𝑌_𝑡 is given by Cov𝑌_𝑡−𝑖,𝑌_𝑡−𝑗=𝜎²_𝜉𝛽₁^{|𝑖−𝑗|}1−𝛽²₁+𝜎_𝜖𝜉𝛽₂1𝛽₁−𝛼₁×𝛽₁^{|𝑖−𝑗|}1−𝛽²₁−𝛼₁^{|𝑖−𝑗|}1−𝛼₁𝛽₁+𝛽₁𝛽₁^{|𝑖−𝑗|}1−𝛽²₁1−𝛼₁𝛽₁+𝜎²_𝜖𝛽²₂𝛽₁−𝛼₁1−𝛼₁𝛽₁×𝛽₁𝛽₁^{|𝑖−𝑗|}1−𝛽²₁−𝛼₁𝛼₁^{|𝑖−𝑗|}1−𝛼²₁,∀_{𝑖,𝑗∈ℤ}. (4)

Considering that 𝛼₁=𝛽₁, we have Cov𝑌_𝑡−𝑖,𝑌_𝑡−𝑗=𝛽₁^{|𝑖−𝑗|}1−𝛽²₁𝜎²_𝜉+2𝜎_𝜖𝜉𝛽₂×𝛽₁1−𝛽²₁+𝑖−𝑗𝛽₁+𝜎²_𝜖𝛽²₂1−𝛽²₁𝑖−𝑗+1+𝛽²₁1−𝛽²₁,for𝑖,𝑗∈ℤ. (5)

In regard to the structure of covariance between the stochastic processes 𝑋_𝑡 and 𝑌_𝑡, for 𝛼₁≠𝛽₁, we have Cov𝑌_𝑡−𝑖,𝑋_𝑡−𝑗=𝜎_𝜖𝜉𝛼₁^{|𝑖−𝑗|}1−𝛼₁𝛽₁+𝜎²_𝜖𝛼₁𝛽₂𝛼₁^{|𝑖−𝑗|}1−𝛼₁𝛽₁1−𝛼²₁,∀_{𝑖,𝑗∈ℤ}. (6)

When 𝛼₁=𝛽₁, the covariance function under study takes the form Cov𝑌_𝑡−𝑖,𝑋_𝑡−𝑗=𝜎²_𝜖𝛽₂𝛼₁𝛼₁^{|𝑖−𝑗|}1−𝛼²₁²+𝜎_𝜖𝜉𝛼₁^{|𝑖−𝑗|}1−𝛼²₁,∀_{𝑖,𝑗∈ℤ}. (7)

By writing out the stochastic system of (1) in matrix notation, the bivariate stochastic process 𝐙_𝑡=[𝑋_𝑡𝑌_𝑡]^′ can be expressed as 𝐙_𝑡=𝛼₀𝛽₀+𝛼₁0𝛽₂𝛽₁𝑋_𝑡−1𝑌_𝑡−1+𝜖_𝑡𝜉_𝑡=𝐜+Φ₁𝐙_𝑡−1+𝝐_𝑡,𝑡∈ℤ, (8) where 𝝐_𝑡=[𝜖_𝑡𝜉_𝑡]^′ is the 2-dimensional Gaussian white noise random vector.

Hence, at each date 𝑡,𝑡∈ℤ, the conditional stochastic process 𝐙_𝑡|_{𝐙𝑡−1=𝐳𝑡−1} follows a bivariate Gaussian distribution, 𝐙_𝑡|_{𝐙𝑡−1=𝐳𝑡−1}⌢𝒩₂(𝝁_{𝑡|𝑡−1},𝛀_{𝑡|𝑡−1}),where the two-dimensional conditional mean value vector and the variance-covariance matrix are, respectively, given by 𝝁_{𝑡|𝑡−1}=𝐜+Φ₁𝐙_𝑡−1,𝛀_{𝑡|𝑡−1}𝜎≡𝛀=^2𝜖𝜎_𝜖𝜉𝜎_𝜖𝜉𝜎²_𝜉. (9)

Straightforward computations lead us to the following factoring of the probability density function of 𝐙_𝑡 conditional to 𝐙_𝑡−1=𝐳_𝑡−1: 𝑓_{𝐙𝑡|𝐙𝑡−1}𝐳_𝑡∣𝐳_𝑡−1=𝑓_{𝑋𝑡|𝐙𝑡−1}𝑥_𝑡∣𝐳_𝑡−1×𝑓_{𝑌𝑡|𝑋𝑡,𝐙𝑡−1}𝑦_𝑡∣𝑥_𝑡,𝐳_𝑡−1. (10)

Thus, the joint distribution of the pair 𝑋_𝑡 and 𝑌_𝑡 conditional to the values of the process at the previous date 𝑡−1, 𝐙_𝑡−1, can be decomposed into the product of the marginal distribution of 𝑋_𝑡|_𝐙𝑡−1 and the conditional distribution of 𝑌_𝑡|_{𝑋𝑡,𝐙𝑡−1}. Both densities follow univariate Gaussian probability laws: 𝑋_𝑡∣_{𝐙𝑡−1=𝐳𝑡−1}⌢𝒩𝛼₀+𝛼₁𝑥_𝑡−1,𝜎²_𝜖,foreachdate𝑡,𝑡∈ℤ. (11) Also, 𝑌_𝑡|_{𝑋𝑡=𝑥𝑡,𝐙𝑡−1=𝐳𝑡−1} follows a Gaussian distribution with 𝐸𝑌_𝑡|_{𝑋𝑡=𝑥𝑡,𝐙𝑡−1=𝐳𝑡−1}=𝛽₀+𝛽₁𝑦_𝑡−1+𝛽₂𝑥_𝑡−1+𝜎_𝜖𝜉𝜎²_𝜖𝑥_𝑡−𝛼₀−𝛼₁𝑥_𝑡−1=𝜓₀+𝜓₁𝑥_𝑡+𝜓₂𝑥_𝑡−1+𝛽₁𝑦_𝑡−1, (12) where 𝜓₁=(𝜎_𝜖𝜉/𝜎²_𝜖) or, for interpretive purposes, 𝜓₁=(𝜎_𝜉/𝜎_𝜖)𝜌_𝜖𝜉. The parameter 𝜓₁ describes, thus, a weighted correlation between the error terms 𝜖_𝑡 and 𝜉_𝑡. The weight corresponds to the ratio of their standard deviations. Moreover, 𝜓₀=𝛽₀−𝜓₁,𝜓₂=𝛽₂−𝜓₁𝛼₁.

The variance has the following structure: Var𝑌_𝑡∣𝑋_𝑡=𝑥_𝑡,𝐙_𝑡−1=𝐳_𝑡−1=𝜎²_𝜉−𝜎²_𝜖𝜉𝜎²_𝜖=𝜎²_𝜉1−𝜌²_𝜖𝜉≡𝜓₃. (13)

The conditional distribution of 𝑌_𝑡|_{𝑋𝑡,𝐙𝑡−1} can be interpreted as a straight-line relationship between 𝑌_𝑡 and 𝑋_𝑡,𝑋_𝑡−1, and 𝑌_𝑡−1. Additionally, it is worth mentioning that if 𝜌_𝜖𝜉=±1 or 𝜎²_𝜉=0, the above conditional distribution degenerates into its mean value. Henceforth, we will discard these particular cases, which means that 𝜓₃≠0.

3. Factoring the Likelihood Based on Monotone Missing Data Pattern

We focus here on theoretical background for factoring the likelihood function from the VAR(1) Model when there are missing values in the data. Suppose that we have the following monotone pattern of missing data: 𝑥₀𝑥₂⋯𝑥_𝑚−1𝑥_𝑚⋯𝑥_𝑛−1𝑦₀𝑦₂⋯𝑦_𝑚−1. (14)

That is, there are 𝑛 observations available from the stochastic process {𝑋_𝑡}_𝑡∈ℤ, whereas due to some uncontrolled factors it was only possible to record 𝑚(𝑚<𝑛) observations from the stochastic process {𝑌_𝑡}_𝑡∈ℤ. In other words, there are 𝑛−𝑚 missing observations from 𝑌_𝑡.

Let the observed bivariate sample of size 𝑛 with missing values: 𝑥₀,𝑦₀,𝑥₁,𝑦₁,…,𝑥_𝑚−1,𝑦_𝑚−1,𝑥_𝑚,…,𝑥_𝑛−1; (15) denote a realisation of the random process 𝐙_𝑡=[𝑋_𝑡𝑌_𝑡]^′,𝑡∈ℤ, which follows a vectorial autoregressive model of order one. The likelihood function, 𝐿(𝜽), is given by 𝐿(𝜽)≡𝑓_{𝐙0,𝐙1,…,𝐙𝑚−1,𝑋𝑚,…,𝑋𝑛−1}×𝐳₀,𝐳₁,…,𝐳_𝑚−1,𝑥_𝑚,…,𝑥_𝑛−1=𝑓_𝐙0𝐳₀^𝑚−1_𝑡=1𝑓_{𝐙𝑡|𝐙𝑡−1}𝐳_𝑡∣𝐳_𝑡−1𝑓_{𝑋𝑚|𝐙𝑚−1}𝑥_𝑚∣𝐳_𝑚−1×^𝑛−1_𝑡=𝑚+1𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡∣𝑥_𝑡−1;𝜽=𝑓_𝐙0𝐳₀^𝑚−1_𝑡=1𝑓_{𝐙𝑡|𝐙𝑡−1}𝐳_𝑡∣𝐳_𝑡−1×^𝑛−1_𝑡=𝑚𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡∣𝑥_𝑡−1, (16) where 𝜽=[𝛼₀𝛼₁𝜎²_𝜖𝛽₀𝛽₁𝛽₂𝜎²_𝜉𝜎_𝜖𝜉]^′ is the 8-dimensional vector of population parameters. To lighten notation, we assume that there is no need for conditioning the arguments of the above probability density functions on the values of the processes at date 𝑡−1. The likelihood function becomes 𝐿(𝜽)=𝑓_𝐙0𝐳₀^𝑚−1_𝑡=1𝑓_{𝐙𝑡|𝐙𝑡−1}𝐳_𝑡^𝑛−1_𝑡=𝑚𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡. (17)

Two points must be emphasised: first, we emphasise that the maximum likelihood estimators (m.l.e.) for the unknown vector of parameters will be obtained by maximising the natural logarithm of the above likelihood function. Second, a worthwhile improvement in reducing the complexity of the function to maximise is to determine the conditional maximum likelihood estimators regarding the first pair of random variables, 𝐙₀=[𝑋₀𝑌₀]^′, as deterministic and maximising the log-likelihood function conditioned on the values 𝑋₀=𝑥₀ and 𝑌₀=𝑦₀. The loss of efficiency of the estimators obtained from such a procedure is negligible when compared with the exact maximum likelihood estimators computed by iterative techniques. Even for moderate sample sizes, the first pair of observations makes a negligible contribution to the total likelihood. Hence, the exact m.l.e. and the conditional m.l.e. turn out to have the same large sample properties, Hamilton [13]. Hereafter, we restrict the study to the conditional loglikelihood function.

Despite the above solutions for reducing the complexity of the problem, some difficulties still remain. The loglikelihood equations are intractable. To go over this problem we have to factorise the conditional likelihood function. From (17) we get 𝐿(𝜽)=^𝑚−1_𝑡=1𝑓_{(𝑋𝑡,𝑌𝑡)|)(𝑋𝑡−1,𝑌𝑡−1}𝑥_𝑡,𝑦_𝑡^𝑛−1_𝑡=𝑚𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡=^𝑚−1_𝑡=1𝑓_{𝑋𝑡|)(𝑋𝑡−1,𝑌𝑡−1}𝑥_𝑡×𝑓_{𝑌𝑡|𝑋𝑡),(𝑋𝑡−1,𝑌𝑡−1}𝑦_𝑡×^𝑛−1_𝑡=𝑚𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡=^𝑛−1_𝑡=1𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡^𝑚−1_𝑡=1𝑓_{𝑌𝑡|𝑋𝑡),(𝑋𝑡−1,𝑌𝑡−1}𝑦_𝑡. (18)

So as to work out the analytical expressions for the unknown parameters under study, we have to decompose the entire likelihood function (18) into easily manipulated components.

For the Gaussian VAR processes, the conditional maximum likelihood estimators coincide with the least squares estimators [13]. Therefore, we may find a solution to the problem just raised in the geometrical context. The identification of such components relies on two of the most famous theorems in the Euclidean space: the Orthogonal Decomposition Theorem and the Approximation Theorem [14, Volume I, pages 572–575]. Based on these tools it is straightforward to establish that the estimation subspaces associated with the conditional distributions 𝑋_𝑡|_𝑋𝑡−1 and 𝑌_𝑡|_{𝑋𝑡,𝑋𝑡−1,𝑌𝑡−1} are, by construction, orthogonal to each other. This means that each element belonging to one of those subspaces is uncorrelated with each element that pertains to their orthogonal complement. Hence, events that happen on one subspace provide no information about events on the other subspace.

The aforementioned arguments guarantee that the decomposition of the joint likelihood in two components can be carried out with no loss of information for the whole estimation procedure. From (18) we can, thus, decompose the conditional loglikelihood function as follows: 𝑙≡𝑙(𝜽)=log𝐿(𝜽)=^𝑛−1_𝑡=1log𝑓_{𝑋𝑡|𝑋𝑡−1}𝑥_𝑡+^𝑚−1_𝑡=1log𝑓_{𝑌𝑡|𝑋𝑡),(𝑋𝑡−1,𝑌𝑡−1}𝑦_𝑡=𝑙₁+𝑙₂. (19)

Henceforth, 𝑙₁ denotes the loglikelihood from the marginal distribution of 𝑋_𝑡, based on the whole sampled data with dimension 𝑛, that is, 𝑥₀,𝑥₁,…,𝑥_𝑛−1. The function 𝑙₂ represents the loglikelihood from the conditional density of 𝑌_𝑡|_{𝑋𝑡,𝐙𝑡−1} computed by the bivariate sample of size 𝑚: 𝑥₀,𝑦₀,𝑥₁,𝑦₁,…,𝑥_𝑚−1,𝑦_𝑚−1. (20)

The components 𝑙₁ and 𝑙₂ of (19) will be maximised separately in Section 4.1.

4. Maximum Likelihood Estimators for the Parameters

In Section 4.1 the m.l.e. of the parameters from the fragmentary VAR(1) Model are deduced. The precision of the estimators is examined in Section 4.2.

4.1. Analytical Expressions

Theoretical developments carried out in this section rely on solving the loglikelihood equations obtained from the factored loglikelihood given by (19). Before proceeding with theoretical matters, we introduce some relevant notation in the ensuing paragraphs.

Let 𝑋_𝑘^(𝑙)=(1/𝑘)∑^𝑘_𝑡=1𝑋_𝑡−𝑙 represent the sample mean lagged 𝑙 time units, 𝑙=0,1. The subscript 𝑘,𝑘=1,…,𝑛−1, allows us to identify the number of observations that takes part in the computation of the sample mean. A similar notation is used for denoting the sample mean of the random sample 𝑌₀,…,𝑌_𝑘, for 𝑘=1,…,𝑚−1, 𝑌_𝑘^(𝑙). According to this new definition, the sample variance of each univariate random variable based on 𝑘 observations and lagged 𝑙 time units is denoted by 𝛾^(𝑙)_𝑋,𝑘=1𝑘^𝑘_𝑡=1𝑋_𝑡−𝑙−𝑋_𝑘^(𝑙)2,𝛾^(𝑙)_𝑌,𝑘=1𝑘^𝑘_𝑡=1𝑌_𝑡−𝑙−𝑌_𝑘^(𝑙)2,𝑙=0,1. (21)

Let 𝛾^∗_𝑋,𝑘(1)=(1/𝑘)∑^𝑘_𝑡=1(𝑋_𝑡−𝑋_𝑘⁽⁰⁾)(𝑋_𝑡−1−𝑋_𝑘⁽¹⁾) describe the sample autocovariance coefficient at lag one for the stochastic process 𝑋_𝑡, based on 𝑘 observations. Its counterpart for the stochastic process 𝑌_𝑡, 𝛾^∗_𝑋,𝑘(1), is obtained by changing notation accordingly. The sample autocorrelation coefficient of the random process 𝑋_𝑡 at lag one is denoted by 𝜌⁽¹⁾_𝑋,𝑘=𝛾^∗_𝑋,𝑘(1)/𝛾⁽⁰⁾_𝑋,𝑘𝛾⁽¹⁾_𝑋,𝑘. The empirical covariance between the random processes 𝑋_𝑡 and 𝑌_𝑡 lagged one time unit is represented by 𝛾^∗_𝑋𝑌(1)=1𝑚−1^𝑚−1_𝑡=1𝑋_𝑡−𝑋⁽⁰⁾_𝑚−1𝑌_𝑡−1−𝑌⁽¹⁾_𝑚−1,𝛾forlaggedvalueson𝑌,_∗𝑌𝑋(1)=1𝑚−1^𝑚−1_𝑡=1𝑋_𝑡−1−𝑋⁽¹⁾_𝑚−1𝑌_𝑡−𝑌⁽⁰⁾_𝑚−1,forlaggedvalueson𝑋. (22)

The sample covariance coefficient of 𝑋_𝑡 and 𝑌_𝑡 computed from 𝑙 time units lag for each series is given by 𝛾^(𝑙)_𝑋𝑌=1𝑚−1^𝑚−1_𝑡=1𝑋_𝑡−𝑙−𝑋^(𝑙)_𝑚−1𝑌_𝑡−𝑙−𝑌^(𝑙)_𝑚−1,with𝑙=0,1. (23)(i)Maximising the loglikelihood function 𝑙₁: Using the results (11) and (19), we readily find the following m.l.e. 𝛼₀=𝑋⁽⁰⁾_𝑛−1−𝛼₁𝑋⁽¹⁾_𝑛−1,𝛼₁=𝛾^∗_{𝑋,𝑛−1}(1)𝛾⁽¹⁾_{𝑋,𝑛−1},𝜎_𝜖²=𝑆𝑆_𝑅,𝑛−1 (24) where 𝑆𝑆_𝑅 is the respective residual sum of squares.(ii)Maximising the loglikelihood function 𝑙₂: Based on (12) and (13) we get the loglikelihood function 𝑙₂𝑙₂=^𝑚−1_𝑡=1log𝑓_{𝑌𝑡|𝑋𝑡,𝑋𝑡−1,𝑌𝑡−1}𝑦_𝑡=−𝑚−12()−log2𝜋𝑚−12log𝜓₃−12𝜓₃×^𝑚−1_𝑡=1𝑦_𝑡−𝜓₀−𝜓₁𝑥_𝑡−𝜓₂𝑥_𝑡−1−𝛽₁𝑦_𝑡−1². (25)

We readily find out that the m.l.e. for the parameters under study are given by 𝜓₀=𝑌⁽⁰⁾_𝑚−1−𝜓₁𝑋⁽⁰⁾_𝑚−1−𝜓₂𝑋⁽¹⁾_𝑚−1−𝛽₁𝑌⁽¹⁾_𝑚−1,𝜓₁=1𝛾⁽⁰⁾_{𝑋,𝑚−1}𝛾⁽⁰⁾_𝑋𝑌−𝜓₂𝛾^∗_{𝑋,𝑚−1}(1)−𝛽₁𝛾^∗_𝑋𝑌(1,𝜓₂)=11−𝜌⁽¹⁾_{𝑋,𝑚−1}²𝛾⁽¹⁾_{𝑋,𝑚−1}×𝛾^∗_𝑌𝑋(1)−𝛾⁽⁰⁾_𝑋𝑌𝛾^∗_{𝑋,𝑚−1}(1)𝛾⁽⁰⁾_{𝑋,𝑚−1}−𝛽₁𝛾⁽¹⁾_𝑋𝑌+𝛽₁𝛾^∗_𝑋𝑌(1)𝛾^∗_{𝑋,𝑚−1}(1)𝛾⁽⁰⁾_{𝑋,𝑚−1},𝜓₃=𝑆𝑆^∗_𝑅,𝛽𝑚−1₁=1𝛾⁽¹⁾_{𝑌,𝑚−1}+𝑌⁽¹⁾_𝑚−1²×𝛾^∗_{𝑌,𝑚−1}(1)+𝑌⁽⁰⁾_𝑚−1𝑌⁽¹⁾_𝑚−1−𝜓₀𝑌⁽¹⁾_𝑚−1−𝜓₁𝛾^∗_𝑋𝑌(1)−𝜓₁𝑋⁽⁰⁾_𝑚−1𝑌⁽¹⁾_𝑚−1−𝜓₂𝛾⁽¹⁾_𝑋𝑌−𝑋⁽¹⁾_𝑚−1𝑌⁽¹⁾_𝑚−1, (26) where 𝑆𝑆^∗_𝑅 denotes the corresponding residual sums of squares.

Using the results from Section 2 we get the following estimators for the original parameters: 𝛽₀=𝜓₀+𝜓₁𝛼₀,𝛽₂=𝜓₂+𝜓₁𝛼₁,𝜎_𝜖𝜉=𝜓₁𝜎²_𝜖,𝜎²_𝜉=𝜓₃+𝜎²_𝜖𝜉𝜎²_𝜖. (27)

Thus, the analytical expressions for the estimators of the mean values, variances, and covariances of the VAR(1) Model are given by 𝜇_𝑋=𝛼₀1−𝛼₁,𝜇_𝑌=𝛼₀𝛽₂+𝛽₀1−𝛼₁1−𝛼₁1−𝛽₁,𝜎²_𝑋=𝜎²_𝜖1−𝛼²₁,𝜎²_𝑌=𝜎²_𝜉1−𝛽²₁+2𝜎_𝜖𝜉𝛽₁𝛽₂1−𝛼₁𝛽₁1−𝛽²₁+𝜎²_𝜖𝛽²₂1+𝛼₁𝛽₁1−𝛼²₁1−𝛽²₁1−𝛼₁𝛽₁𝛼₁≠𝛽₁,𝜎²_𝑌=𝜎²_𝜉1−𝛼²₁𝜎+2^𝜖𝜉𝛼₁𝛽₂1−𝛼²₁²+𝜎²_𝜖𝛽²₂1+𝛼²₁1−𝛼²₁³,𝜎_𝑋𝑌=𝜎_𝜖𝜉1−𝛼₁𝛽₁+𝛼₁𝛽₂𝜎²_𝜖1−𝛼₁𝛽₁1−𝛼²₁(atthesamedate𝑡,𝑡∈ℤ). (28)

These estimators will play a central role in the following sections.

4.2. Precision of the Estimators

In the section, the precision of the maximum likelihood estimators underlying equations (28) is derived. The whole analysis will be separated in three stages. First, we study the statistical properties of the vector Θ, where Θ=[Θ₁Θ₂]^′, with Θ₁=[𝛼₀𝛼₁𝜎²_𝜖]^′ and Θ₂=[𝜓₀𝜓₁𝜓₂𝛽₁𝜓₃]^′. For notation consistency, the unknown parameter 𝛽₁ is either denoted by 𝛽₁ or 𝜓₄. That is, 𝜓₄≡𝛽₁. Secondly, we derive the precision of the m.l.e. of the original parameters of the VAR(1) Model (see (1)). Finally, we will focus our attention on the estimators for the mean vector and the variance-covariance matrix at lag zero of the VAR(1) model with a monotone pattern of missingness.

There are a few points worth mentioning. From Section 3 we know that there is no loss of information in maximising separately the loglikelihood functions 𝑙₁ and 𝑙₂ (19). As a consequence, the variance-covariance matrix associated with the whole set of estimated parameters is a block diagonal matrix. For sufficiently large sample size, the distribution of the maximum likelihood estimator is accurately approximated by the following multivariate Gaussian distribution: Θ≈𝒩₈Θ₁Θ₂,𝐈₁⁻¹𝟎𝟎𝐈₂⁻¹, (29) where 𝐈₁ and 𝐈₂ denote the Fisher information matrices, respectively, from the components 𝑙₁ and 𝑙₂ of the loglikelihood function (see (19)). There is an asymptotic equivalence between the Fisher information matrix and the Hessian matrix (see [8, ch.2]). Moreover, as long as Θ→Θ there is also an asymptotic equivalence between the Hessian matrix computed at the points Θ and Θ. Henceforth, the Fisher information matrices from (29) are estimated, respectively, by 𝐈₁̂=−𝜕²𝑙₁𝜕Θ₁𝜕Θ^′_1Θ1=Θ₁,𝐈₂̂𝜕=−²𝑙₁𝜕Θ₂𝜕Θ^′_2Θ2=Θ₂. (30)

To lighten notation, from now on we suppress the “hat” from the consistent estimators of the information matrices.

The variance-covariance matrix for Θ₁ takes the following form: 𝐈₁⁻¹=𝜎²_𝜖(𝑛−1)𝛾⁽¹⁾_{𝑋,𝑛−1}(0)×𝛾⁽¹⁾_{𝑋,𝑛−1}(0)+𝑋⁽¹⁾_𝑛−1²−𝑋⁽¹⁾_𝑛−10−𝑋⁽¹⁾_𝑛−110002𝜎²_𝜖𝛾⁽¹⁾_{𝑋,𝑛−1}(0). (31)

We stress that there is orthogonality between the error and the estimation subspaces underlying the loglikelihood function 𝑙₁.

Calculating the second derivatives of the loglikelihood function 𝑙₂ results in the following approximate information matrix:𝐈₂=1𝜓₃𝑚−1(𝑚−1)𝑋_𝑚−1(𝑚−1)𝑋⁽¹⁾_𝑚−1(𝑚−1)𝑌⁽¹⁾_𝑚−10(𝑚−1)𝑋∑_𝑚−1^𝑚−1_𝑡=1𝑋∑²_𝑖^𝑚−1_𝑡=1𝑋_𝑡𝑋∑_𝑡−1^𝑚−1_𝑡=1𝑋_𝑡𝑌_𝑡−10(𝑚−1)𝑋∑⁽¹⁾_𝑚−1^𝑚−1_𝑡=1𝑋_𝑖𝑋∑_𝑖−1^𝑚−1_𝑡=1𝑋∑²_𝑡−1^𝑚−1_𝑡=1𝑋_𝑡−1𝑌_𝑡−10(𝑚−1)𝑌∑⁽¹⁾_𝑚−1^𝑚−1_𝑡=1𝑋_𝑖𝑌∑_𝑖−1^𝑚−1_𝑡=1X_𝑡−1𝑌∑_𝑡−1^𝑚−1_𝑡=1𝑌²_𝑡−100000𝑚−12𝜓₃. (32)

Once again, we mention that there is orthogonality between the error and the estimation subspaces underlying the loglikelihood function 𝑙₂. The matrix 𝐈₂ can be written in a compact form: 𝐈₂=1𝜓₃𝐈₂₁𝟎𝟎𝐼₂₂, (33) where the (4×4) submatrix 𝐈₂₁ and the scalar 𝐼₂₂ are, respectively, defined as 𝐈₂₁=𝐔^′𝐔,𝐼₂₂=𝑚−12𝜓₃², (34) with 𝐔=1𝑋₁𝑋₀𝑌₀1𝑋₂𝑋₁𝑌₁⋮⋮⋮⋮1𝑋_𝑚−1𝑋_𝑚−2𝑌_𝑚−2. (35)

Using the above partition of 𝐈₂ it is rather simple to compute the inverse matrix. In fact, 𝐈₂⁻¹=𝜓₃𝐈^−𝟏₂₁𝟎𝟎𝐼⁻¹₂₂, (36) with 𝐈⁻¹₂₁=(𝐔^′𝐔)⁻¹ and 𝐼⁻¹₂₂=(2/(𝑚−1))(𝜓₃)².

Unfortunately, there is no explicit expression for the inverse matrix 𝐈⁻¹₂₁. As a result, there are no explicit expressions for the approximate variance-covariance of the m.l.e. for the vector of unknown parameters Θ₂.

Now, we have to analyse the precision of the m.l.e. of the original parameters of the VAR(1) Model, that is, 𝚼=[𝛼₀𝛼₁𝜎²_𝜖𝛽₀𝛽₁𝛽₂𝜎²_𝜉𝜎_𝜖𝜉]^′.

Recalling from Section 2, the one-one monotone functions that relate the vector of parameters under consideration, that is, Θ₂=𝜓₀𝜓₁𝜓₂𝜓₄𝜓₃^′,𝚼₂=𝛽₀𝛽₁𝛽₂𝜎²_𝜉𝜎_𝜖𝜉^′𝜓,are₀=𝛽₀−𝛼₀𝜓₁,𝜓₁=𝜎_𝜖𝜉𝜎²_𝜖,𝜓₂=𝛽₂−𝜓₁𝛼₁𝜓₃=𝜎²_𝜉−𝜎²_𝜖𝜓²₁,𝜓₄≡𝛽₁. (37)

The parameters 𝛼₀,𝛼₁, and 𝜎²_𝜖 remain unchanged. A key assumption in the following developments is that neither the estimates of the unknown parameters nor the true values fall on the boundary of the allowable parameter space.

The variance-covariance matrix of the m.l.e. for the vector of parameters 𝚼 is obtained by the first-order Taylor expansion at 𝚼. We also use the chain rule for derivatives of vector fields ([for details, see [14, Volume II, pages 269–275]).

Writing the vector of parameters 𝚼 as a function of the vector Θ, the respective first-order partial derivatives can be joined together in the following partitioned matrix: 𝐃𝐃=¹𝐃₂𝐃₃𝐃₄, (38) where the (3×3) submatrix 𝐃₁ corresponds to the first-order partial derivatives of the vector 𝚼₁≡Θ₁=[𝛼₀𝛼₁𝜎²_𝜖]^′ with respect to itself, which means that 𝐃₁ is nothing but the identity matrix of order 3, 𝐃₁=𝐈₃. On the other hand, this statement also means that the derivatives of the parameters under consideration with respect to either 𝜓₀,𝜓₁,𝜓₂,𝜓₃, or 𝜓₄ are zero. In other words, the (3×5) submatrix 𝐃₂ is equal to the null vector, that is, 𝐃₂=𝟎.

The (5×3) submatrix 𝐃₃ and the (5×5) submatrix 𝐃₄ are composed by the first-order partial derivatives of each component of the vector of parameters 𝚼₂=[𝛽₀𝛽₁𝛽₂𝜎²_𝜉𝜎_𝜖𝜉]^′ with respect to, respectively, 𝛼₀,𝛼₁,𝜎²_𝜖 and 𝜓₀,𝜓₁,𝜓₂,𝜓₄,𝜓₃. Their structures are, thus, given by 𝐃₃=𝜓₁000000𝜓₁000𝜓²₁00𝜓₁,𝐃₄=1𝛼₀000000100𝛼₁10002𝜓₁𝜎²_𝜖0010𝜎²_𝜖000=𝐃₄₁𝐃₄₂𝐃₄₃𝐷₄₄. (39)

For finding out the approximate variance-covariance matrix of the maximum likelihood estimators for the unknown vector of parameters 𝚼, it is only necessary to pre- and postmultiply the variance-covariance matrix arising from expressions (29), (31), and (36) by, respectively, the matrix 𝐃 and its transpose, 𝐃^′. More precisely, 𝚺_𝚼≈𝐃𝐈⁻¹𝐃^′=𝐈₃𝟎_3×5𝐃₃𝐃₄𝐈₁⁻¹𝟎𝟎𝐈₂⁻¹𝐈₃𝐃^′₃𝟎_5×3𝐃^′₄. (40) Hence, 𝚺_𝚼≈𝐈₁⁻¹𝐈₁⁻¹𝐃^′₃𝐈₁⁻¹𝐃^′₃^′𝐃₃𝐈₁⁻¹𝐃^′₃+𝐃₄𝐈₂⁻¹𝐃^′₄, (41) with 𝚺_𝚼 denoting the variance-covariance matrix of the m.l.e. for the vector of unknown parameters 𝚼. A more detailed analysis of the variance-covariance matrix (41) can be found in Nunes [12, ch.3, p.91-92].

We can now deduce the approximate variance-covariance matrix of the maximum likelihood estimators for the mean vector and the variance-covariance matrix at lag zero of the VAR(1) Model with a monotone pattern of missingness, represented by Ξ=[𝛼₀𝛼₁𝜎²_𝜖𝜇_𝑋𝜇_𝑌𝜎²_𝑋𝜎²_𝑌𝜎_𝑋𝑌]^′. The first-order partial derivatives of the vector Ξ with respect to the vector 𝚼 are placed in a matrix that is denoted by 𝐅. It takes the following form: 𝐅𝐅=¹𝐅₂𝐅₃𝐅₄. (42)

According to the partition of the matrix 𝐃 into four blocks—expression (38)—we partition the matrix 𝐅 into the following blocks: the (3×3) submatrix 𝐅₁ corresponds to the partial derivatives of 𝛼₀,𝛼₁, and 𝜎²_𝜖 with respect to themselves. As a consequence, 𝐅₁ is the identity matrix of order 3, that is, 𝐅₁=𝐈₃. Regards to the (3×5) sub-matrix 𝐅₂, its elements correspond to the partial derivatives of 𝛼₀,𝛼₁, and 𝜎²_𝜖 with respect to 𝛽₀,𝛽₁,𝛽₂,𝜎²_𝜉, and 𝜎_𝜖𝜉. Therefore, 𝐅₂=𝟎. The partial derivatives of 𝜇_𝑋,𝜇_𝑌,𝜎²_𝑋,𝜎²_𝑌, and 𝜎_𝑋𝑌 with respect to 𝛼₀,𝛼₁, and 𝜎²_𝜖 are gathered together in the (5×3) sub-matrix 𝐅₃: 𝐅₃=𝑓³₁₁𝑓³₁₂0𝑓³₂₁𝑓³₂₂00𝑓³₃₂𝑓³₃₃0𝑓³₄₂𝑓³₄₃0𝑓³₅₂𝑓³₅₃, (43) where 𝑓³₁₁=11−𝛼₁,𝑓³₁₂=𝛼₀1−𝛼₁²,𝑓³₂₁=𝛽₂1−𝛽₁1−𝛼₁,𝑓³₂₂=𝛼₀𝛽₂1−𝛽₁1−𝛼₁²,𝑓³₃₂=2𝛼₁𝜎²_𝜖1−𝛼²₁²,𝑓³₃₃=11−𝛼²₁,𝑓³₄₂=2𝛽₂𝜎_𝜖𝜉𝛽²₁1−𝛼²₁²+𝜎²_𝜖𝛽₂𝛽₁1−𝛼²₁+𝛼₁1−𝛼²₁𝛽²₁×1−𝛼₁𝛽₁1−𝛼²₁²1−𝛽²₁⁻¹,𝑓³₄₃=𝛽²₂1+𝛼₁𝛽₁1−𝛼²₁1−𝛽²₁1−𝛼₁𝛽₁,𝑓³₅₂=𝜎_𝜖𝜉𝛽₁1−𝛼₁𝛽₁²+𝜎²_𝜖𝛽₂1+𝛼²₁1−2𝛼₁𝛽₁1−𝛼₁𝛽₁²1−𝛼²₁²,𝑓³₅₃=𝛼₁𝛽₂1−𝛼₁𝛽₁1−𝛼²₁. (44)

The 5-dimensional square sub-matrix 𝐅₄ corresponds to the partial derivatives of 𝜇_𝑋,𝜇_𝑌,𝜎²_𝑋,𝜎²_𝑌, and 𝜎_𝑋𝑌 with respect to 𝛽₀,𝛽₁,𝛽₂,𝜎²_𝜉,𝜎_𝜖𝜉: 𝐅₄=𝑓00000₄₂₁𝑓⁴₂₂𝑓⁴₂₃00000000𝑓⁴₄₂𝑓⁴₄₃𝑓⁴₄₄𝑓⁴₄₅0𝑓⁴₅₂𝑓⁴₅₃0𝑓⁴₅₅, (45) with its nonnull elements taking the following analytical expressions: 𝑓⁴₂₁=11−𝛽₁,𝑓⁴₂₂=𝛽₀1−𝛽₁²+𝛼₀𝛽₂1−𝛼₁1−𝛽₁²,𝑓⁴₂₃=𝛼₀1−𝛼₁1−𝛽₁,𝑓⁴₄₂=2𝜎²_𝜉𝛽₁1−𝛽²₁²+2𝜎_𝜖𝜉𝛽₂1+𝛽²₁1−2𝛼₁𝛽₁1−𝛼₁𝛽₁²1−𝛽²₁²+2𝜎²_𝜖𝛽²₂𝛼₁1−𝛽²₁+𝛽₁1−𝛼²₁𝛽²₁1−𝛼²₁1−𝛽²₁²1−𝛼₁𝛽₁²,𝑓⁴₄₃=21−𝛼₁𝛽₁1−𝛽²₁𝛽₁𝜎_𝜖𝜉+𝜎²_𝜖𝛽₂1+𝛼₁𝛽₁1−𝛼²₁,𝑓⁴₄₄=11−𝛽²₁,𝑓⁴₄₅=2𝛽₁𝛽₂1−𝛼₁𝛽₁1−𝛽²₁,𝑓⁴₅₂=1−𝛼²₁𝜎_𝜖𝜉𝛼₁+𝜎²_𝜖𝛼²₁𝛽₂1−𝛼²₁1−𝛼₁𝛽₁²,𝑓⁴₅₃=𝜎²_𝜖𝛼₁1−𝛼₁𝛽₁1−𝛼²₁,𝑓⁴₅₅=11−𝛼₁𝛽₁. (46)

Straightforward calculations have paved the way to the desired partitioned variance-covariance matrix, called here 𝚺_Ξ, 𝚺_Ξ≈𝐅𝚺_𝚼𝐅^′≈𝐅𝐃𝐈⁻¹𝐃^′𝐅^′=𝚺_Ξ¹¹𝚺_Ξ¹²𝚺_Ξ²¹𝚺_Ξ²², (47) with its submatrices defined by 𝚺_Ξ¹¹=𝐈₁⁻¹,𝚺_Ξ¹²=𝐈₁⁻¹𝐅₃+𝐅₄𝐃₃^′=𝐈₁⁻¹𝐅^′₃+𝐃^′₃𝐅^′₄,𝚺_Ξ²¹=𝐅₃+𝐅₄𝐃₃𝐈₁⁻¹=𝚺_𝐹^12′,𝚺_Ξ²²=𝐅₃𝐈₁⁻¹𝐅^′₃+𝐅₄𝐃₃𝐈₁⁻¹𝐅^′₃+𝐅₃𝐈₁⁻¹𝐃^′₃𝐅^′₄+𝐅₄𝐃₃𝐈₁⁻¹𝐃^′₃+𝐃₄𝐈₂⁻¹𝐃^′₄𝐅^′₄=𝐅₃+𝐅₄𝐃₃𝐈₁⁻¹𝐅₃+𝐅₄𝐃₃^′+𝐅₄𝐃₄𝐈₂⁻¹𝐅₄𝐃₄^′=𝐆𝐈₁⁻¹𝐆^′+𝐇𝐈₂⁻¹𝐇^′. (48)

The matrix 𝐆 that has just been defined as 𝐆=𝐅₃+𝐅₄𝐃₃ corresponds to the first-order partial derivatives from the composite functions that relate 𝜇_𝑋,𝜇_𝑌,𝜎²_𝑋,𝜎²_𝑌, and 𝜎_𝑋𝑌 with the vector of parameters Θ₁. The elements of the matrix 𝐇=𝐅₄𝐃₄ are the first-order partial derivatives from the composite functions that relate 𝜇_𝑋,𝜇_𝑌,𝜎²_𝑋,𝜎²_𝑌, and 𝜎_𝑋𝑌 with the vector of unknown parameters Θ₂.

The 3-dimensional square sub-matrix 𝚺_Ξ¹¹ corresponds to the approximate covariance structure between the m.l.e. of the parameters 𝛼₀,𝛼₁, and 𝜎²_𝜖. The (3×5) sub-matrix 𝚺_Ξ¹² is composed of the approximate covariances between the m.l.e. that have just been cited and 𝜇_𝑋,𝜇_𝑌,𝜎²_𝑋,𝜎²_𝑌, and 𝜎_𝑋𝑌; its transpose is denoted by 𝚺_Ξ²¹. This is the reason why 𝚺_Ξ¹², or 𝚺_Ξ²¹, results from the product of the variance-covariance matrix 𝐈₁⁻¹ and 𝐆. The 5-dimensional square sub-matrix 𝚺_Ξ²² is formed by the covariances between the m.l.e. for 𝜇_𝑋,𝜇_𝑌,𝜎²_𝑋,𝜎²_𝑌, and 𝜎_𝑋𝑌.

The main point of the section is to study the variances and covariances that take part of the sub-matrix 𝚺_Ξ²². Thus, it is of interest to further explore its analytical expression. The matrix 𝐆 takes a cumbersome form. The most efficient way to deal with it is to consider its partition rather than the whole matrix at once.

Let 𝐆=𝐅₃+𝐅₄𝐃₃=𝐆₁₁𝐆₁₂𝐆₂₁𝐺₂₂, (49) where the (4×2) sub-matrix 𝐆₁₁ takes the form 𝐆₁₁=11−𝛼₁𝛼₀1−𝛼₁²𝜓₁1−𝛽₁+𝛽₂1−𝛼₁1−𝛽₁𝛼₀1−𝛼₁1−𝛽₁𝜓₁+𝛽₂1−𝛼₁02𝛼₁𝜎²_𝜖1−𝛼²₁²0𝑔¹¹₄₂, (50)

with 𝑔¹¹₄₂=𝜓₁𝑓⁴₄₃+2𝛽₂1−𝛽²₁1−𝛼₁𝛽₁²1−𝛼²₁²×𝜎_𝜖𝜉𝛽²₁1−𝛼²₁²+𝜎²_𝜖𝛽₂𝛽₁1−𝛼²₁+𝛼₁1−𝛼²₁𝛽²₁, (51) where 𝑓⁴₄₃ is defined by (45).

The 4-dimensional column vector 𝐆₁₂, the 2-dimensional row vector 𝐆₂₁ and the scalar 𝐺₂₂ are, respectively, given by 𝐆₁₂=1001−𝛼²₁𝜓²₁1−𝛽²₁+𝛽₂1−𝛼₁𝛽₁1−𝛽²₁×2𝜓₁𝛽₁+𝛽₂1+𝛼₁𝛽₁1−𝛼²₁^′,𝐆₂₁=011−𝛼₁𝛽₁×𝜎_𝜖𝜉𝛽₁1−𝛼₁𝛽₁+𝜎²_𝜖1−𝛼²₁×𝜓₁𝛼₁+𝛽₂1+𝛼²₁1−2𝛼₁𝛽₁1−𝛼²₁1−𝛼₁𝛽₁,𝐺₂₂=11−𝛼₁𝛽₁𝜓₁+𝛼₁𝛽₂1−𝛼²₁. (52)

On the other hand, we can also make the following partition of the matrix 𝐇=𝐅₄𝐃₄: 𝐇𝐇=¹¹𝐇₁₂𝐇₂₁𝐻₂₂, (53) where the sub-matrix 𝐇₁₁ corresponds to the first order partial derivatives of the vector [𝜇_𝑋𝜇_𝑌𝜎²_𝑋𝜎²_𝑌]^′ with respect to the vector [𝜓₀𝜓₁𝜓₂𝜓₄]^′, whereas their derivatives with respect to the parameter 𝜓₃ constitute the sub-matrix 𝐇₁₂. The sub-matrix 𝐇₂₁ is composed of the first order partial derivatives of 𝜎_𝑋𝑌 with respect to each component of the vector [𝜓₀𝜓₁𝜓₂𝜓₄]^′. Finally, the scalar 𝐻₂₂=𝜕𝜎_𝑋𝑌/𝜕𝜓₃=0.

The desired variance-covariance matrix can therefore be written in the following partitioned form: 𝚺_Ξ²²=𝚺₁²²𝚺₂²²𝚺₃²²𝚺₄²², (54) with 𝚺₁²²=𝜎²_𝜖𝐆₁₁𝐔^′_𝑅𝐔_𝑅⁻¹𝐆^′₁₁+2𝜎⁴_𝜖𝐆𝑛−1¹²𝐆^′₁₂+𝜓₁𝐇₁₁𝐔^′𝐔⁻¹𝐇^′₁₁+2𝜓²₃𝐇𝑚−1¹²𝐇^′₁₂,𝚺₂²²=𝜎²_𝜖𝐆₁₁𝐔^′_𝑅𝐔_𝑅⁻¹𝐆^′₂₁+2𝜎⁴_𝜖𝐆𝑛−1¹²𝐆^′₂₂+𝜓₃𝐇₁₁𝐔^′𝐔⁻¹𝐇^′₂₁,𝚺₃²²=𝜎²_𝜖𝐆₂₁𝐔^′_𝑅𝐔_𝑅⁻¹𝐆^′₁₁+2𝜎⁴_𝜖𝐆𝑛−1²²𝐆^′₁₂+𝜓₃𝐇₂₁𝐔^′𝐔⁻¹𝐇^′₁₁,𝚺₄²²=𝜎²_𝜖𝐆₂₁𝐔^′_𝑅𝐔_𝑅⁻¹𝐆^′₂₁+2𝜎⁴_𝜖𝐆𝑛−1²²𝐆^′₂₂+𝜓₃𝐇₂₁𝐔^′𝐔⁻¹𝐇^′₂₁, (55) where the matrix 𝐔 is defined by (35). The matrix 𝐔_𝑅 takes the form 𝐔_𝑅=1𝑋₀1𝑋₁⋮⋮1𝑋_𝑚−2. (56)

In short, the matrix defined by (54) corresponds to the approximate variance-covariance matrix of the m.l.e. for the mean vector and variance-covariance matrix at lag zero for the VAR(1) Model with missing data. We cannot write down explicit expressions for those variances and covariances. The limitation arises from the inability to invert the matrix product 𝐔^′𝐔 in analytical terms (see (36)). Hence, its inverse can only be accomplished by numerical techniques using the observed sampled data. This point will be pursued further in Section 5.

Despite the above restrictions, several investigations can be done regarding the amount of additional information obtained by making full use of the fragmentary data available. The strength of the correlation between the stochastic processes here plays a crucial role. These ideas will be developed in Section 5.

5. Simulation Studies

In this section, we analyse the effects of using different strategies to estimate the mean value of the stochastic process {𝑌_𝑡,𝑡∈ℤ}, denoted by 𝜇_𝑌. More precisely, the bivariate modelling scheme and its univariate counterparts are compared. Simulation studies are carried out to evaluate the relative efficiency of the estimators with interest.

The m.l.e. of the mean value of the stochastic process {𝑌_𝑡,𝑡∈ℤ} based on the VAR(1) Model is obtained by the second equation of the system of (28). We need to compare this estimator to those obtained by considering the univariate stochastic process {𝑌_𝑡,𝑡∈ℤ} itself. More precisely, having in mind that we are handling a bivariate VAR(1) Model, the corresponding marginal model is the ARMA(2,1) [10, 11]. On the other hand, the AR(1) Model is one of the most popular models due to its practical importance in time series modelling. Therefore, the behaviour of the AR(1) Model will be also evaluated. In short, we will compare the performance of the VAR(1) Model with both the ARMA(2,1) and the AR(1) Models.

To avoid any confusion between the parameters coming from the bivariate and the univariate modelling strategies, from now on we denote the parameter from the VAR(1) Model by 𝜇_VAR, whereas those from the ARMA(2,1) and the AR(1) Models are represented by 𝜇_ARMA and 𝜇_AR, respectively.

The bivariate VAR(1) Model is described by the system of (1). Thus, the univariate stochastic process {𝑌_𝑡,𝑡∈ℤ} follows an ARMA(2,1) Model, and the m.l.e. of the mean value are given by 𝜇_ARMA=𝛽₀1−𝛼₁1−𝛽₁−𝛽₁. (57)

On the other hand, if we assumed that {𝑌_𝑡,𝑡∈ℤ} followed an AR(1) Model, the m.l.e. of the mean value would be given by 𝜇_AR=𝛽₀1−𝛽₁. (58)

Next, we will compare the performance of the estimators (57) and (58) with the m.l.e. based on the VAR(1) Model (second equation of the system (28)). It is important to stress that the strategy behind the AR(1) Model has not taken into account the relationship between the stochastic processes {𝑋_𝑡,𝑡∈ℤ} and {𝑌_𝑡,𝑡∈ℤ}. This feature will certainly introduce an additional noise in the overall estimation procedure.

Following the techniques used in Section 4.2 for determining the precision of the estimators under consideration, here we also have used the first-order Taylor expansion at the mean value 𝜇_𝑌 for computing the estimate of the variance of 𝜇_𝑌.

Considering the ARMA(2,1) Model, let 𝜽=[𝛽₀𝛽₁𝛼₁]^′ be the vector of the unknown parameters. Then, Var𝜇_ARMA≈³_𝑖=1𝜕𝜇_ARMA𝜕𝜃_𝑖∣_𝜃𝑖=𝜃_𝑖²Var𝜃_𝑖+2³³_{𝑖=1𝑗=𝑖+1}𝜕𝜇_ARMA𝜕𝜃_{𝑖𝜃𝑖=}𝜃_𝑖×𝜕𝜇_ARMA𝜕𝜃_{𝑗𝜃𝑗=}𝜃_𝑗Cov𝜃_𝑖,𝜃_𝑗. (59)

In regard to the AR(1) Model, 𝜇_AR is given by (58) and Var𝜇_AR≈2𝛽₀1−𝛽₁³Cov𝛽₀,𝛽₁+Var𝛽₀1−𝛽₁²+𝛽²₀Var𝛽₁1−𝛽₁⁴. (60)

Improvements in choosing the sophisticated m.l.e. for 𝜇_𝑌 based on the VAR(1) Model rather than considering its univariate counterparts are next discussed. Simulation studies are carried out to evaluate the relative efficiency of the estimators under consideration.

The data were generated by the VAR(1) Model (system of (1)). In order to make comparisons on the same basis, a few assumptions to the parameters of the VAR(1) Model are made. We consider that 𝜇_𝑋=𝜇_𝑌=0. These restrictions have no influence on the results because they are equivalent to 𝛼₀=𝛽₀=0, that is, the constant terms of the VAR(1) Model are equal to zero (system of (1)). Additionally, we introduce the restriction 𝜎²_𝜖=𝜎²_𝜉=1.

Since the correlation coefficient regulates the supply of information between the stochastic processes {𝑋_𝑡}_𝑡∈ℤ and {𝑌_𝑡}_𝑡∈ℤ, particular emphasis is given to this parameter. Using the grid of points 𝜌_𝜖𝜉=0.1,0.5,0.75,0.9, the Gain index is computed. We stress that the value 𝜌_𝜖𝜉=1 is not allowable in this context (see Section 2 for the details).

We analyse the performance of the estimators based on different sample sizes, 𝑛=50,100,250, and 500. The simulations reported next are based on different percentages of missing observations referred to the dimension of the sampled data from the auxiliary random process {𝑋_𝑡}_𝑡∈ℤ. Simulation runs for each combination of the parameters are based on 1000 replicates.

It is worth emphasising that the estimates of the covariance terms that take part of the variances given by (59) and (60) were computed by the R package tseries [15].

The simulation goes as follows: after each simulation run, the relative efficiency of 𝜇_VAR with respect to each estimator 𝜇_ARMA and 𝜇_AR is quantified by the Gain index, GI₁ and GI₂, respectively, expressed as percentage: GI₁=Var𝜇_ARMA−Var𝜇_VARVar𝜇_ARMA×100%,GI₂=Var𝜇_AR−Var𝜇_VARVar𝜇_AR×100%. (61)

A word of notation: the above quantities, that is, GI₁ and GI₂, were computed from the estimates of the corresponding variances. To lighten the notation, we skipped the conventional nomenclature used to represent the estimates.

If GI₁>0, then 𝜇_VAR is more precise than 𝜇_ARMA. Otherwise, 𝜇_VAR loses precision, and 𝜇_ARMA becomes a better estimator for the mean value of {𝑌_𝑡}_𝑡∈ℤ. A similar reasoning applies to the comparison between 𝜇_VAR and 𝜇_AR.

Figures 1 and 2 display the main results from the simulation studies. The estimators 𝜇_VAR and 𝜇_ARMA are compared in Figure 1, whereas Figure 2 exhibits the comparison between 𝜇_VAR and 𝜇_AR. For each combination of the parameters of the model, we represent graphically the gain indexes as functions of the percentage of missing data in the sampled data from the stochastic process {𝑌_𝑡}_𝑡∈ℤ.

(a) 𝜌𝜖𝜉=0.1

(b) 𝜌𝜖𝜉=0.5

(c) 𝜌𝜖𝜉=0.75

(d) 𝜌𝜖𝜉=0.9

Figure 1 

Graphical representation of GI1. The data were obtained from a VAR(1) Model, with 𝛼0=𝛽0=0, 𝛼1=0.6, 𝛽1=0.7, and 𝛽2=0.8.

(a)  𝜌𝜖𝜉=0.1

(b) 𝜌𝜖𝜉=0.5

(c) 𝜌𝜖𝜉=0.75

(d) 𝜌𝜖𝜉=0.9

Figure 2 

Graphical representation of GI2. The data were obtained from a VAR(1) Model, with 𝛼0=𝛽0=0,𝛼1=0.6,𝛽1=0.7, and 𝛽2=0.8.

Either Figures 1 or 2 shows that the plot of the gain index against the percentage of missing data in the sample from the stochastic process {𝑌_𝑡}_𝑡∈ℤ behaves roughly as a linear function, regardless of the combination of the parameters. In outline, the more the percentage of missing values in the sampled data is, the more precise is the estimator 𝜇_VAR when compared with the univariate context, that is, 𝜇_ARMA or 𝜇_AR (see Figures 1 and 2).

Further, the gain in precision by using the sophisticated estimator 𝜇_VAR rather than 𝜇_ARMA or 𝜇_AR increases as the strength of the linear relationship between the processes {𝑋_𝑡}_𝑡∈ℤ and {𝑌_𝑡}_𝑡∈ℤ (described by the correlation coefficient) rises from 𝜌=0.1 to 𝜌=0.9. This statement is true for both the ARMA(2,1) and AR(1) modelling schemes (see Figures 1 and 2).

A final point to highlight from the comparison between Figures 1 and 2 is that the increase in precision obtained by using the estimator for the mean value of {𝑌_𝑡}_𝑡∈ℤ based on the VAR(1) Model is higher when we compare its performance with the results from the AR(1) Model than when we compare the VAR(1) Model with the ARMA(2,1) Model. This feature emphasises the idea that has already been raised that the ARMA(2,1) Model describes more accurately the dynamics of the stochastic process {𝑌_𝑡}_𝑡∈ℤ than the AR(1) Model does. In short, it seems that the AR(1) Model is not a good approach in this context because it incorporates a noise term related to the simulation scheme that we cannot control.

Summing up, the estimator 𝜇_VAR is preferable to those explored in the univariate context, that is, either 𝜇_ARMA or 𝜇_AR.

6. Conclusions

This article deals with the problem of missing data in an univariate sample. We have considered an auxiliary complete data set, whose underlying stochastic process is serially correlated with the former by the VAR(1) Model structure. We have proposed maximum likelihood estimators for the relevant parameters of the model based on a monotone missing data pattern. The precision of the estimators has also been derived. Special attention has been given to the estimator for the mean value of the stochastic process whose sampled data has missing values, 𝜇_𝑌.

We have compared the performance of the estimator for 𝜇_𝑌 based on the VAR(1) Model with a monotone pattern of missing data with those obtained from both the ARMA(2,1) Model and the AR(1) Model. By simulation studies, we have showed that the estimator derived in this article based on the VAR(1) Model performs better than those derived from the univariate context. It is essential to emphasise that, even numerically, it was quite difficult to compute the precision of the later estimators as we have shown in Section 4.2.

A compelling question remains unresolved. From an applied point of view, it would be extremely useful to develop estimators for the dynamics of the stochastic processes. More precisely, we would like to get estimators for the correlation and cross-correlation matrices as well as their precision when there are missing observations in one of the data sets. It was not possible to achieve this goal based on maximum likelihood principles. As we have shown in Section 4.2, we have only developed estimators for the covariance matrix at lag zero. In future research, we will try to solve this problem in the framework of Kalman filter.

Acknowledgments

This work was financed by the Portuguese Foundation for Science and Technology (FCT), Projecto Estratégico PEst-OE/MAT/UI0209/2011. The authors are also thankful for the comments of the two anonymous referees.