Predicting seasonal and hydro-meteorological impact in environmental variables modelling via Kalman filtering

Abstract

This study focuses on the potential improvement of environmental variables modelling by using linear state-space models, as an improvement of the linear regression model, and by incorporating a constructed hydro-meteorological covariate. The Kalman filter predictors allow to obtain accurate predictions of calibration factors for both seasonal and hydro-meteorological components. This methodology can be used to analyze the water quality behaviour by minimizing the effect of the hydrological conditions. This idea is illustrated based on a rather extended data set relative to the River Ave basin (Portugal) that consists mainly of monthly measurements of dissolved oxygen concentration in a network of water quality monitoring sites. The hydro-meteorological factor is constructed for each monitoring site based on monthly precipitation estimates obtained by means of a rain gauge network associated with stochastic interpolation (kriging). A linear state-space model is fitted for each homogeneous group (obtained by clustering techniques) of water monitoring sites. The adjustment of linear state-space models is performed by using distribution-free estimators developed in a separate section.

Appendix

1.1 Distribution-free estimators for the mean and for the transition matrix

In the parameters estimation of state-space models were performed distribution-free estimators developed from the original work by Costa and Alpuim (2010). However, in that work it was proposed a distribution-free estimator for state-space models with univariate observations. Thus, a straightforward generalization of these estimators is presented in order to allow their application to a class of multivariate state-space models that largely covers the present work’s needs.

To estimate the unknown parameters in the model

$$ {\mathbf Y}_t={\mathbf H}_t {\mathbf \beta}_t+{\mathbf e}_t $$

(11)

$$ {\mathbf \beta}_t={\mathbf \mu}+{\mathbf \Upphi}\left({\mathbf \beta}_{t-1}-{\mathbf \mu}\right)+{\mathbf \epsilon}_t $$

(12)

it is assumed a set of observations ${{\mathcal Y}_n=(\mathbf{Y}_1, {\mathbf Y}_2, \ldots, {\mathbf Y}_n),}$ and regular matrices of known constants ${\mathbf H}_1, {\mathbf H}_2, \ldots, {\mathbf H}_n$ are available. The mean vector ${\mathbf \mu}$ can be easily estimated by the method of moments, i.e., $\widehat{\mathbf \mu}=n^{-1}\sum_{t=1}^n {\mathbf H}_t^{-1}{\mathbf Y}_t.$

As variables Y _t are not stationary, we are not under the usual conditions of the consistency of generalized method of moments. Thus, it is necessary to establish additional conditions to guarantee this consistency. By construction, the estimator $\widehat{\mathbf \mu}$ of the mean vector is unbiased, so we can guarantee its consistency by proving that $var(\widehat{\mathbf \mu})\rightarrow {\mathbf 0}$ when $n\rightarrow+\infty,$ and thus establishing sufficient conditions. Covariance matrix of $\widehat{\mathbf \mu}$ is given by

$$ \begin{aligned} var(\widehat{\mathbf \mu}) =\frac{1}{n^2}\sum_{t=1}^n\sum_{s=1}^nE \left[\left(\bf{H}_t^{-1}\bf{Y}_t-{\mathbf \mu}\right)\left({\mathbf H}_s^{-1} \bf{Y}_s-{\mathbf{\mu}}\right)^{\prime}\right]\\ &=\frac{1}{n^2}\sum_{t=1}^n\sum_{s=1}^nE\left[\left({\mathbf \beta}_t-{\mathbf \mu}\right)\left({\mathbf \beta}_s-{\mathbf{\mu}}\right)^{\prime} \right]+\frac{1}{n^2}\sum_{t=1}^n\bf{H}_t^{-1}{\mathbf \Upsigma}_ {\mathbf e} {\mathbf H}_t^{\prime -1}. \end{aligned} $$

Applying the Kronecker product ⊗ and the operator vec, we get

$$ vec\left[var\left(\widehat{{\mathbf{\mu}}}\right)\right]=\frac{1}{n^2}\left [\sum_{t=1}^n\sum_{s=1}^n {\mathbf{\Upgamma}}_{{\mathbf{\beta}}}(|t-s|)\right]+\frac{1} {n^2}\left[\sum_{t=1}^n\left(\bf{H}_t^{-1}\otimes\bf{H}_t^{-1}\right) \right]vec\left({\mathbf{\Upsigma}}_{\mathbf{e}} \right). $$

Under the stationarity conditions of process {β _t } the first parcel is an O _p , seeing that $\sum_{k=-\infty}^{+\infty}\mathbf{\Upgamma}(k)<\infty,$ (e.g., Hamilton 1994, p. 279). To guarantee that

$$ \frac{1}{n^2}\left[\sum_{t=1}^n\left(\bf{H}_t^{-1}\otimes\bf{H}_t^{-1} \right)\right]vec\left({\mathbf{\Upsigma}}_{\mathbf{e}}\right)\buildrel {n\longrightarrow+\infty} \over {\longrightarrow}{\mathbf{0}} $$

it is sufficient to admit the additional condition |h ⁻¹ _t,(i,j) | < c for all $t=1,2,.., i,j=1,2,\ldots,m$ and for some positive constant c, where h ⁻¹ _t,(i,j) represents the (i, j) element of H ⁻¹ _t matrix.

The autoregressive matrix $\mathbf{\Upphi}$ is estimated by means of covariance structure of process ${\mathbf H}_t^{-1}{\mathbf Y}_t.$ We see that

$$ \begin{aligned} {\mathbf{\Upgamma}}_{\mathbf{{H}}^{-1}{\mathbf{Y}}}(k)&=E\left[\left({\mathbf{H}}_ {t+k}^{-1}\bf{Y}_{t+k}-{\mathbf{\mu}}\right)\left(\bf{H}_{t}^{-1}\bf{Y}_{t}- {\mathbf{\mu}}\right)^{\prime}\right]\\ &={\mathbf{\Upgamma}}_{{\mathbf{\beta}}}(k)={\mathbf{\Upgamma}}_k. \end{aligned} $$

In a VAR(1) process, the relation $\mathbf{\Upgamma}_k=\mathbf{\Upphi \Upgamma}_{k-1}$ is valid, for k = 1, 2, .... Thus, we proposed the autoregressive matrix estimator $\widehat{\Upphi}$ based on the least squares method of these equations by taking $k=1,2,\ldots,\ell_\mathbf{\Upphi}.$ Thus, we have

$$ \widehat{{\mathbf{\Upphi}}}=\left(\sum_{k=1}^{\ell_{{\mathbf{\Upphi}}}}\widehat {{\mathbf{\Upgamma}}}_{k+1}\widehat{{\mathbf{\Upgamma}}}^{\prime}_{k} \right)\left(\sum_{k=1}^{\ell_{{\mathbf{\Upphi}}}}\widehat{{\mathbf{\Upgamma}}} _{k}\widehat{{\mathbf{\Upgamma}}}^{\prime}_{k}\right)^{-1} $$

where $\widehat{\mathbf{\Upgamma}}_k=\frac{1}{n}\sum_{t=1}^{n-k}\left[\left(\mathbf {H}_{t+k}^{-1}\mathbf{Y}_{t+k}-\widehat{\mathbf{\mu}}\right)\left(\mathbf {H}_{t}^{-1}\bf{Y}_{t}-\widehat{\mathbf{\mu}}\right)^{\prime}\right].$

By construction, the autoregressive matrix estimator is consistent, since $\widehat{\mathbf{\Upgamma}}_k$ is a consistent estimator of $\mathbf{\Upgamma}_k.$ Whereas we have proposed a consistent estimator to $\mathbf{\mu},$ we consider that the mean vector $\mathbf{\mu}$ is known. To analyse the consistency of $\widehat{\mathbf{\Upgamma}}_k$ we have

$$ \begin{aligned} \widehat{\mathbf \Upgamma}_k&=\frac{1}{n}\sum_{k=1}^{n-k}\left[({\mathbf \beta}_{t+k}-{\mathbf \mu}+{\mathbf H}^{-1}_{t+k}{\mathbf e}_{t+k})({\mathbf \beta}_ {t}-{\mathbf \mu}+\bf{H}^{-1}_{t}\bf{e}_{t})^{\prime}\right]\\ &=\frac{1}{n}\sum_{k=1}^{n-k}\left[({\mathbf \beta}_{t+k}-{\mathbf \mu}) ({\mathbf \beta}_{t}-{\mathbf \mu})^{\prime}+({\mathbf \beta}_{t+k}-{\mathbf \mu}){\mathbf e}^{\prime}_t {\mathbf H^{\prime}} _t^{-1}\right. \\ &\left. +{\mathbf e}^{\prime}_{t+k}{\mathbf H^{\prime}}_{t+k}^{-1}({\mathbf \beta} _{t}-{\mathbf \mu})^{\prime}+{\mathbf e}^{\prime}_{t+k} {{\mathbf{H^{\prime}}}}_{t+k} ^{-1}{\mathbf e}^{\prime}_t{\mathbf H^{\prime}}_t^{-1}\right]. \end{aligned} $$

Under the previously established condition, the last three parcels converge in probability to a null matrix. Indeed, by defining the second parcel as ${\mathbf A}=[A_{ij}]_{i,j=1,2,\ldots,m}$ and, with some algebraic manipulation, we have

$$ A_{ij}=\frac{1}{n}\sum_{t=1}^{n-k}\left(({\mathbf \beta}_{t,(i)}-\mu_i) \sum_{s=1}^me_{t,(s)}h^{-1}_{t,(s,j)}\right) $$

and considering σ _e,(r,s) = cov(e _t,(r), e _t,(r)), the variance is given by

$$ var(A_{ij})=\sigma^2_{{\mathbf \beta}_{i}}\sum_{r=1}^m\sum_{s=1}^m \sigma_{e,(r,s)}\left(\frac{1}{n^2}\sum_{t=1}^{n-k}h^{-1}_{t,(r,j)}h^{-1} _{t,(s,j)}\right) $$

If the additional condition |h ⁻¹ _t,(i,j) | < c is valid, this parcel tends to 0 when $n\rightarrow+\infty.$ In a similar way, we defined the third parcel by $\mathbf{B}=[B_{ij}]_{i,j=1,2,\ldots,m}$ with elements given by

$$ B_{ij}=\frac{1}{n}\sum_{t=1}^{n-k}\left((\beta_{t,(j)}-\mu_j) \sum_{s=1}^me_{t+k,(s)}h^{-1}_{t+k,(i,s)}\right) $$

with variance

$$ var(B_{ij})=\sigma^2_{\beta_j}\sum_{r=1}^m\sum_{s=1}^m\sigma_{e,(r,s)} \left(\frac{1}{n^2}\sum_{t=1}^{n-k}h^{-1}_{t+k,(i,s)}h^{-1}_{t+k,(i,r)} \right) $$

Again, we guarantee that B _ij  = O _p through the same condition |h ⁻¹ _t,(i,j) | < c. As we shall see, this condition is a sufficient condition, as the last parcel also tends to a null matrix. Indeed, if we denote the last parcel as $\mathbf{C}=[C_{ij}]_{i,j=1,2,\ldots,m},$ we have

$$ C_{ij}=e_{t,(j)}\sum_{r=1}^me_{t+k,(r)}h^{-1}_{t+k,(i,r)}h^{-1}_{t,(i,r)} $$

with variance given by

$$ var(C_{ij})=\sigma^2_{e_j}\frac{1}{n^2}\sum_{t=1}^{n-k}\sum_{r=1}^m\sum_ {s=1}^mh^{-1}_{t+k,(i,r)}h^{-1}_{t,(i,r)}h^{-1}_{t+k,(i,s)}h^{-1}_ {t,(i,s)}\sigma_{e,(r,s)}. $$

These results allow us to conclude that if |h ⁻¹ _t,(i,j) | < c, the estimator $\widehat{\mathbf{\Upgamma}}_k$ is consistent to $\mathbf{\Upgamma},$ when we replace the mean vector $\mathbf{\mu}$ by a consistent estimator.

1.2 Distribution-free estimators to noise variances

The estimation of covariance matrices of errors terms e _t and ε _t is an important and difficult step at the same time. At times, the recursive procedures applied to the obtained Gaussian likelihood estimates diverge or produce non-positive semidefined matrices. Sometimes, these problems occur when the initial solution is not as close to estimates as necessary. We propose an estimator to $\mathbf{\Upsigma}_\mathbf{\varepsilon}$ based on covariance structure of a VAR(1) stationary process.

We know that the relation $\mathbf{\Upsigma}_{\mathbf{\beta}}=\mathbf{\Upphi\Upsigma\Upphi}^{\prime}+ \mathbf{\Upsigma}_{\mathbf{\varepsilon}}$ is valid in a VAR(1) stationary process, or by applying the Kronecker product ⊗ and the operator vec

$$ vec\left({\mathbf{\Upsigma}}_{{\mathbf{\beta}}}\right)=\left[{\mathbf{ I}}_{m^2}-({\mathbf{\Upphi}}\otimes {\mathbf{\Upphi}})\right]^{-1}vec\left({\mathbf{\Upsigma}}_ {{\mathbf{\varepsilon}}}\right). $$

By applying the vec operator to the equation $\mathbf{\Upgamma}_k=\mathbf{\Upphi\Upgamma}_{k-1},$ with $k=1,2,\ldots,$ we have:

$$ vec\left({\mathbf \Upgamma}_k\right)=vec\left({\mathbf \Upphi}^k {\mathbf \Upsigma}_{\mathbf \beta}\right) =\left({\mathbf I}_m\otimes {\mathbf \Upphi}^k \right)vec\left({\mathbf \Upsigma}_{\mathbf \beta}\right) =\left({\mathbf I}_m\otimes {\mathbf \Upphi}^k \right)\left[{\mathbf I}_{m^2}-({\mathbf \Upphi}\otimes {\mathbf \Upphi})\right]^{-1}vec\left({\mathbf \Upsigma}_ {\mathbf \varepsilon}\right) $$

or

$$ vec\left({\mathbf \Upgamma}_k\right)^{\prime}=vec\left({\mathbf \Upsigma}_ {\mathbf \varepsilon}\right)^{\prime}\left[{\mathbf I}_{m^2}- ({\mathbf \Upphi}\otimes {\mathbf \Upphi})^{\prime}\right]^{-1}\left({\mathbf I}_m\otimes {\mathbf \Upphi}^k \right)^{\prime}. $$

Note that the matrix ${\mathbf \Upsigma}_{\mathbf \varepsilon}$ is symmetric, that is, $vec\left({\mathbf \Upsigma}_{\mathbf \varepsilon}\right)^{\prime}_{1,(j-1) m+i-1}=vec \left({\mathbf \Upsigma}_{\mathbf \varepsilon}\right)^{\prime}_{1,im+j}$ with 1 ≤ i ≤ m − 1 and 1 ≤ j ≤ i. Thus, we constructed a line matrix $vec\left({\mathbf \Upsigma}_{\mathbf \varepsilon}\right)^{\prime*},$ with m + m(m − 1)/2 columns, that we got from $vec\left({\mathbf \Upsigma}_{\mathbf \varepsilon}\right)^{\prime}$ by removing the elements $vec\left({\mathbf \Upsigma}_{\mathbf \varepsilon}\right)^{\prime}_{1,im+j},$ with 1 ≤ i ≤ m − 1 and 1 ≤ j ≤ i.

By applying the same methodology to the matrix ${\mathbf \Updelta}_k$ defined as

$$ {\mathbf \Updelta}_k=\left[{\mathbf I}_{m^2}-({\mathbf \Upphi}\otimes {\mathbf \Upphi})^{\prime}\right]^{-1}\left({\mathbf I}_m\otimes {\mathbf \Upphi}^k \right)^{\prime}, $$

we summed the columns (two by two) with the index im + j and (j − 1)m + i − 1, with 1 ≤ i ≤ m − 1 and 1 ≤ j ≤ i, thus obtaining a new matrix ${\mathbf \Updelta}_k^*$ with m + m(m − 1)/2 columns.

The estimator for ${\mathbf \Upsigma}_{\mathbf \varepsilon}$ is constructed via the least squares method applied to equations

$$ vec\left({\mathbf \Upgamma}_k\right)^{\prime}=vec\left({\mathbf \Upsigma}_ {\mathbf \varepsilon}\right)^{\prime*} {\mathbf \Updelta}_k^* $$

with $k=1,2,\ldots,\ell_{\mathbf \varepsilon}.$ Thus, we obtained the estimator

$$ vec\left(\widehat{\mathbf \Upsigma}_{\mathbf \varepsilon}\right)^{\prime*} =\left(\sum_{k=1}^{\ell_{\mathbf \varepsilon}}\widehat{\mathbf \Upgamma} _k\widehat{\mathbf \Updelta}_k^{\prime*}\right)\left(\sum_{k=1}^ {\ell_{\mathbf \varepsilon}}\widehat{\mathbf \Updelta}_k^*\widehat {\mathbf \Updelta}_k^{\prime*}\right)^{-1}. $$

(13)

The consistency of $\widehat{\mathbf \Upsigma}_{\mathbf \varepsilon}$ is guaranteed under the same conditions of the consistency of $\widehat{\mathbf \Updelta}_k.$ As we have seen, a sufficient condition for this is |h ⁻¹ _t,(i,j) | < c.

In order to estimate the covariance matrix ${\mathbf \Upsigma}_{\mathbf e},$ we defined

$$ {\mathbf \Uppsi}=\frac{1}{n}\sum_{t=1}^n\left[\left({\mathbf H}_t^{-1} {\mathbf Y}_t-{\mathbf \mu}\right)\left({\mathbf H}_t^{-1}{\mathbf Y}_t-{\mathbf \mu} \right)^{\prime} \right]. $$

Therefore, we had the expectation

$$ \begin{aligned} E({\mathbf \Uppsi})&=\frac{1}{n}\sum_{t=1}^nE\left[\left({\mathbf \beta}_t- {\mathbf \mu}+{\mathbf H}_t^{-1}{\mathbf e}_t\right)\left({\mathbf \beta}_t- {\mathbf \mu}+{\mathbf H}_t^{-1}{\mathbf e}_t\right)^{\prime}\right]\\ &={\mathbf \Upsigma}_{\mathbf \beta}+\frac{1}{n}\sum_{t=1}^n\left({\mathbf H}_t^{-1}{\mathbf \Upsigma}_{\mathbf e}{\mathbf H}_t^{\prime-1}\right). \end{aligned} $$

By applying the vec operator, and with some algebraic manipulation, we got

$$ vec\left(\widehat{\mathbf \Upsigma}_{\mathbf e}\right)^{\prime}=n\left [vec\left(\widehat{\mathbf \Uppsi}\right)^{\prime}-vec\left(\widehat{\mathbf \Upsigma}_{\mathbf \beta}\right)^{\prime}\right]\left[\sum_{t=1}^n\left ({\mathbf H}_t^{-1}\otimes {\mathbf H}_t^{-1}\right)^{\prime}\right]^{-1}. $$

As the matrix ${\mathbf \Upsigma}_{\mathbf e}$ is symmetric, it is necessary to adopt the same procedure as in the estimation of ${\mathbf \Upsigma}_{\mathbf \varepsilon}.$ Thus, we estimated the m + m(m − 1)/2 elements of the covariance matrix.

If we have a consistent estimator to ${\mathbf \Upsigma}_{\mathbf \beta},$ for example given by the proposed estimators to ${\mathbf \Upphi}$ and ${\mathbf \Upsigma}_{\mathbf \varepsilon},$ the consistency of $\widehat{\mathbf \Upsigma}_{\mathbf e}$ boils down to the limit of variance of each element of $vec({\mathbf \Upupsilon})= n vec(\widehat{\mathbf \Uppsi})^{\prime}[\sum_{t=1}^n({\mathbf H}_t^{-1}\otimes {\mathbf H}_t^{-1})^{\prime}]^{-1}.$ The variance of the (i, j) element of ${\mathbf \Upupsilon}$ is given by

$$ n^2a_{ij}^2var\left[\frac{1}{n}\sum_{t=1}^n\left(\beta_{t,i}-\mu_i- \sum_{k=1}^mh^{-1}_{t,(i,k)}e_{t,k}\right)\left(\beta_{t,j}-\mu_j- \sum_{k=1}^mh^{-1}_{t,(j,k)}e_{t,k}\right)\right] $$

where h ⁻¹ _t,(i,j) represents the (i,j) element of the matrix ${\mathbf H}_t^{-1}$ and a _ij the (i, j) element of the matrix $[\sum_{t=1}^n({\mathbf H}_t^{-1}\otimes {\mathbf H}_t^{-1})^{\prime}]^{-1}.$

For simplicity, we adopt β _t,i  − μ _i  = β ^* _t,i . If we take in account that the states ${\mathbf \beta}_t$ are uncorrelated to noise ${\mathbf e}_s$ for all t and s, the previous expression can be decomposed into four parcels. The first parcel has the form

$$ \begin{aligned} &a_{ij}^2var\left(\sum^n_{t=1}\beta_{t,i}^*\beta_{t,j}^*\right)+ a_{ij}^2var\left(\sum_{t=1}^n\beta_{t,i}^*\sum_{s=1}^mh^{-1}_ {t,(j,s)}e_{t,s}\right)\\ &+a_{ij}^2var\left(\sum_{t=1}^n\beta_{t,j}^*\sum_{k=1}^mh^{-1}_ {t,(i,k)}e_{t,k}\right)+a_{ij}^2var\left[\sum_{t=1}^n\left(\sum_{k=1} ^mh^{-1}_{t,(i,k)}e_{t,k}\sum_{s=1}^mh^{-1}_{t,(j,k)}e_{t,s}\right)\right] \end{aligned} $$

The first parcel can be decomposed into

$$ a_{ij}^2var\left(\sum^n_{t=1}\beta_{t,i}^*\beta_{t,j}^*\right)= a_{ij}^2\sum_{t=1}^nvar \left(\beta^*_{t,i}\beta^*_{t,j}\right)+\sum_{t=1}^n\sum_ {^{s=1}_{s \neq t}}^n cov \left(\beta^*_{t,i}\beta^*_{t,j},\beta^*_{s,i}\beta^*_ {s,j}\right) $$

but we can write

$$ var \left(\beta_{t,i}^*\beta_{t,j}^*\right)=cov \left(\beta_{t,i}^{2},\beta_{t,j}^{2}\right)+\sigma^2 _{\beta_i}\sigma^2_{\beta_j}-\gamma^2_{ij}. $$

In order for this parcel to be an O _p , it is sufficient to admit the additional regularity conditions, such as cov(β _t,i β _t,j, β _s,i β _t,j ) for all t and s, that do not depend on time.

The cross terms have a similar structure. For example, the first term can be computed by,

$$ \begin{aligned} a_{ij}^2var\left(\sum_{t=1}^n\beta_{t,i}^*\sum_{s=1}^mh^{-1}_{t,(j,s)}e_{t,s}\right) &=a_{ij}^2\sum_{t=1}^n var\left(\sum_{k=1}^mh^{-1}_{t,(i,k)}\beta^*_{t,j}e_{t,k}\right)\\ &=a_{ij}^2\sum_{t=1}^n\sum_{k=1}^m h^{-2}_{t,(i,k)}\sigma^2_{\beta_j}\sigma^2_{e_k}\\ &=a_{ij}^2\sigma^2_{\beta_j}\sum_{k=1}^m\sigma^2_{e_k}\sum_{t=1}^nh^{-2}_{t,(i,k)}.\\ \end{aligned} $$

So, if we admit that the elements of matrix $\mathbf{H}^{-1}_t$ are limited as c ₁ < |h ⁻¹ _t,(i,j) | < c ₂, where c ₁ and c ₂ are positive constants, it follows that this term is an O _p . In addition to these conditions on h ⁻¹ _t,(i,j) , if we ensure that the vector of error e _t is stationary of fourth-order, then we conclude that the last parcel of variance of the (i, j) element of ${\mathbf \Upupsilon}$ is an O _p , too.

Thus, under the additional stationarity conditions of fourth-order on the vector of disturbances and the above restrictions on the elements of the matrices ${\mathbf H}_t^{-1},$ the proposed distribution-free estimator to ${\mathbf \Upsigma}_{\mathbf e}$ is consistent.