0
$\begingroup$

I want to analyze non-stationary time-series, with the aim of doing a Granger causality analysis. For this i need reliable and time-varying VAR model coefficients and I think using a Kalman filter formalism, I should be able to obtain time varying coefficients. I'm basing my model on an article by Oya, et al, 2007: https://www.ncbi.nlm.nih.gov/pubmed/17184906, and am writing in R.

My issue is with understanding how to obtain an a priori covariance matrix for the estimation error $\Sigma_{t,t-1}$, and how to define the process noize Q. Furthermore, when i run some time-series through my written algorithm, I feel like it is way too accurate, almost like over-fitting. Is this possible and why?

The model is:

$X(t) = F(t)X(t-1)+W(t)$

$Y(t) = H(t) X(t) + E(t)$

Here $F(t)$ is the transition matrix taken to be a diagonal matrix, meaning the coefficients evolve as a random-walk. $H(t)$ contains the lagged values of observations $H(t) = I_m \otimes [Y_{t-1},...,Y_{t-p}]$, with $m$ being the dimensions (2 in my case, since i look at two time-series) and $p$ is the model order. The VAR coefficients are stored in $X(t)$ as $X(t)=Vec[A_1(t),...,A_p(t)]$. Here in the article there appears to be a mistake as the Vec operator is written Vech, which does not reduce to the corresponding VAR model for two time-series of order $p$.

In the article, the process noise covariance matrix is determined as $E[W_t W_t^T]=Q(t)$ and observation noise is $E[E_t E_t^T]=R(t)$. I determine $E(t)$ for the calculation of $R(t)$ as the residuals $E(t) = Y(t) - H(t)\hat{X}(t)$, where $\hat{X}(t)$ is the estimated VAR coefficients. $Q(t)$ i just take as a unit matrix times a small value, $Q(t) = 10^{-3} I$.

My code:

library(fUnitRoots)
library(urca)
library(vars)
library(aod)
library(zoo)
library(tseries)
library(aTSA)

#Im running a test on data available in R as default, the beavers data set

b1<-beaver1
b2<-beaver2

#Since the dimensions don't conform, ill just take 100 observations of both beavers body temperatures.  

sL <- 100
Y<-matrix(0,m,sL)
Y[1,]<-b2$temp[1:sL]
Y[2,]<-b1$temp[1:sL]

#I select the the model order to be *p=3* and define the dimensions as *m=2*, as i have two time-series. 

m<-2
p<-3

#In the cited article, the matrix F is a ((m^2)*p) diagonal matrix. I understand this to define 
#that the coefficients evolve in time as a random-walk. 

F <- diag((m^2)*p) 

#Now I need the a-priori coefficients, or state-estimates for the Kalman filter. As i have no
#idea how to estimate them, i just make a suitable matriks of ones, saying that the first 
#prediction for my state is just the sum of lagged state observations. I call it "X-hat-minus"

Xhatminus <- matrix(1,(m^2)*p,1)

#Now I define the error terms of the coefficients *W* and var process *E*. 
#In Oya, *et al*, 2007, there seems to be a mistake in their dimensions. They must be vectors of 
# length ((m^2)*p), not matrixes of size ((m^2)*p). 
#I choose *W* to be a sequence of random variables in the beginning. 
#Q I understand to be the coefficient estimate error covariance matrix. From the article, i
#understood it to be a unit matrix with all elements 10^-3.5  

Q <- diag((m^2)*p)*10^(-3.5)
E<-matrix(rnorm((m), mean = 0, sd = 1),nrow=(m))
#From E, i understand that the covariance matrix is found as:
R<-E%*%t(E)

#Now i define the *a priori* covariance matrix of the estimation error, calling it *Pmin*, as
#the sigma symbol is hard to write in R. As i have no idea how to find it, i make it a 
#diagonal matrix of appropriate size. 

Pmin = diag((m^2)*p) 

#I store the calculated coefficients for later use in a matrix "X-hat-storage".

Xhatstorage <- matrix(0,(m^2)*p,sL)

#I also store the predicted time-series values in a "Y-hat" matrix, for comparison later with 
#actual data

Yhat<-matrix(0,m,(sL-p))

#I define channels (since im basing it off an EEG analysis article) for convenience in the loop,
# which are essentially the time-series indicators 

ch1<-1
ch2<-2

#Now what follows is my understanding of the Kalman filter, applied to a VAR model, with the 
#observations being the time-series values and the estimation being the coefficients. 

for(j in (p+1):(sL-p)){

#I create the *H* matrix, containing lagged measurements of Y.    

H <- diag(m) %x% t(vec(t(Y[c(ch1,ch2),(j-p):(j-1)])))

#Now i calculate the Kalman gain matrix

K <- (Pmin %*% t(H)) %*% solve((H%*%Pmin%*%t(H) + R)) 

#Here i find the updated coefficients, based on the *a-priori* coefficient estimates  

Xhatplus <- F%*%(Xhatminus + K%*%(Y[c(ch1,ch2),j]-H%*%Xhatminus))

#I now update the coefficient covariance matrix *P*. This corresponds to, in many sources, 
#P_(t|t), calculated from P_(t|t-1)

P <- Pmin - K%*%H%*%Pmin

#And now i calculate the *a-priori* coeficient covariance matrix estimate for the next loop.

Pplus <- (F%*% P %*% F) + Q

#I update the coefficients to use them in the next loop as a-priori. 

Xhatminus <- Xhatplus

#I update the error terms with more properly estimated coefficients. 

E <- Y[c(ch1,ch2),j]-H%*%Xhatminus
R <- E%*%t(E)

#I store the calculated coefficients outside the loop for later use. 

Xhatstorage[,j]<- Xhatplus 

#Also, i store a predicted variable matrix to compare in a plot later with the original data

Yhat[,j]<-H%*%Xhatplus

#I set the new coefficient covariance matrix to be *a-priori* for the next loop. 

Pmin <- Pplus
}

Now I have the coefficients, but they seem to be way too accurate, as seen in the plot below, produced as follows:

plot(1:(sL-p),Y[1,1:(sL-p)])
lines(1:(sL-p),Yhat[1,1:(sL-p)],col="red")

I plot the original data for the first beaver and then the one obtained by the multiplication of the lagged data values and coefficients given by the Kalman filter (inside Yhat).

Points represent the temperature data from beaver1 and red line represents model estimates

I hope this example is quickly and easily reproducible in any R version. I'll reiterate my questions:

1) How do i define the a priori covariance matrix for the estimation error, noted as $\Sigma_{t,t-1}$ in the article and "Pmin" in my code.

2) How do i define the state-noise covariance matrix Q

3) Is the model over-fitting? The estimated values and data line up too well. In fact, the errors are in magnitudes of 10^-19 in some cases.

$\endgroup$
1
$\begingroup$

1) The a priori covariance matrix of the vector state can be dictated by subject matter knowledge, as can the a priori mean vector. Quite often you do not have this detailed subject matter knowledge, in which case you set the initial mean vector to zeroes and the a priori covariance matrix to a diagonal with very large variances (to reflect prior ignorance). A more sophisticated way to proceed is a so-called diffuse prior: some packages in R (like KFAS) make provision for it.

2) If you have no subject matter knowledge to set $Q$ you can always estimate it.

3) If your algorithm tracks the observations it is probably because $Q$ contains variances too large as compared to those in $H$. Variances in $H$ are for the observation error, those in $Q$ are variances in the random process which drives the state vector. If $H << Q$ what you are basically saying is that your observation error is far smaller than the uncertainty introduced by the process error in going from $t$ to $(t+1)$ and the Kalman filter will track almost exactly your trusted observations.

This is a large subject for a short answer such as this. I would recommend you to look at books such as Durbin-Koopman or Shumway-Stoffer (here is the official free pdf of the newest edition) for details.

$\endgroup$
8
  • 1
    $\begingroup$ I would also recommend to use one of the packages in R for state-space estimation (like dlm, KFAS, or MARSS, the last being close to problems such as the one you present). Kalman filtering can be tricky numerically, and these packages may save you a headache or two. $\endgroup$ – F. Tusell Apr 20 '20 at 15:25
  • $\begingroup$ Thank you for the input! In regards to the 3rd point - if I give the Q too large values and this results in the Kalman filter predicting more accurately.. Is there a simple way to put it for a non-econometrics schooled person, why this is bad? What am i missing with this type of over-fitting, when the values for t are always predicted from the lagged values.. doesn't this mean my coefficients are very good? $\endgroup$ – Castle Apr 20 '20 at 15:30
  • 1
    $\begingroup$ Your model overfits when it "learns" your sample and adapts to it. What you want (I guess) is your model to generalize to observations that were not seen in the training sample, only then you can be confident that the model captures the underlying phenomenon. In general, in non-engineering problems where $Q$ might be known, I think you should estimate it. $\endgroup$ – F. Tusell Apr 20 '20 at 15:39
  • $\begingroup$ What i want the most is to have a reliable set of regression coefficients across the time-points available (discounting the first var order lenght of points). These coefficients then i can use to uncover the contributions from input variables over time. In this case, wouldnt saying "my observation error is far smaller than the uncertainty introduced by the process error, going from t to (t+1)" be a good thing? As the Kalman Filter tracks the known observations strongly, giving me in the form of the coefficients or the states accurate information about contributing lags in the variables? $\endgroup$ – Castle Apr 20 '20 at 15:48
  • $\begingroup$ No, it would not necessarily be a good thing. You can always track exactly your data, you only have to assume a random walk hypothesis and near zero observation error! But you probably want to learn about some underlying mechanism driving your data. You should set $H$ according to whatever the true error of observation is. Pretending your observations are more precise than they really are, will make your coefficients fit the signal... and the noise. $\endgroup$ – F. Tusell Apr 20 '20 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.