0
$\begingroup$

I'm posting this again, since after my last post i've been able to advance the code quite alot. I'm still trying to write some code in R to reproduce the model i found in this article.

The idea is to model the signal as a VAR model, but fit the coefficients by a Kalman-filter model. As i understand this would enable me to deal with non-stationary data to a degree.

The model to track the coefficients is:

X(t) = F(t) X(t− 1) +W(t)

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

where H(t) is the Kronecker product between lagged measurements in my time-series Y and a unit vector, and X(t) fills the role of regression-coefficients. F(t) is taken to be an identity matrix, as that should mean we assume coefficients to evolve as a random walk.

In the article, from W(T), the state noise covariance matrix Q(t) is chosen at 10^-3 at first and then fitted based on some iteration scheme. From E(t) the state noise covariance matrix is R(t) substituted by the covariance of the noise term unexplained by the model: Y(t) - H(t)Xhat(t)

I have the a priori covariance matrix of estimation error (denoted Σ in the article) written as P (based on other sources) and the a posteriori as Pmin, since it will be used in the next recursion as a priori, if that makes sense.

So far i've written the following, based on the articles Appendix A 1.2

Y <-                                  *my timeseries, for test purposes two channels of 3000 points*
F <- diag(8)                          # F is (m^2*p by m^2 *p) where m=2 dimensions and p =2 lags
H <- diag(2) %x% t(vec(Y[,1:2]))      #the Kronecker product of vectorized lags Y-1 and Y-2
Xhatminus <- matrix(1,8,1)            # an arbitrary *a priori* coefficient matrix
Q <- diag(8)%x%(10**-7)               #just a number a really low number diagonal matrix, found it used in some examples
R<- 1                                 #Didnt know what else to put here just yet
Pmin = diag(8)                        #*a priori* error estimate, just some 1-s...

Now should start the reccursion. To test i just took the first 3000 points of one trial of my data.

   Xhatstorage <- matrix(0,8,3000)

   for(j in 3:3000){
            H <- diag(2) %x% t(vec(Y[,(j-2):(j-1)]))
            K <- (Pmin %*% t(H)) %*% solve((H%*%Pmin%*%t(H) + R)) ##Solve gives inverse matrix ()^-1
            P <- Pmin - K%*% H %*% Pmin
            Xhatplus <- F%*%( Xhatminus + K%*%(Y[,j]-H%*%Xhatminus) )
            Pplus <- (F%*% P %*% F)  + Q
            Xhatminus <- Xhatplus
            Xhatstorage[,j] <- Xhatplus 
            Pmin <- Pplus
    }

I extracted Xhatplus values into a storage matrix and used them to write this primitive VAR model with them:

Yhat<-array(0,3000)
for(t in 3:3000){
Yhat[t]<- t(vec(Y[,(t-2)])) %*% Xhatstorage[c(1,3),t] + t(vec(Y[,(t-1)])) %*% Xhatstorage[c(2,4),t]
}

The plot looks like

test plot of kalman-filtered coefficient VAR model and actual data

The blue line is VAR with Kalman filter found coefficients, Black is original data.. Only showing the first variable plot here, the other one fits a bit better actually.

I'm having issue understanding how i can better evaluate my coefficients? Why is it so off?

Besides this, how should i better choose the first a priori and a posteriori estimates to start the recursion? Most places i pieced this together from start from arbitrary 0 assumptions in toy models.

Lastly, is this recursion even a correct implementation of Oya et al describe in the article?

$\endgroup$
  • $\begingroup$ Hi: I can't comment on the specific code for the VAR application but your question seems more general in the sense that you are asking:: Are there are there methodologies for estimatiing the system and observation variances in both the A) "classical" version of the KF and the B) "bayesian" version. For A) check out the prediction error decomposition approach of Harvey in his blue book titled: "structural models and the KF". ( his 1990 book ). For B), check out the yellow book by West and Harrison. I think the title of the latter it dynamic linear models. . $\endgroup$ – mark leeds Apr 7 at 15:55

Your Answer

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

Browse other questions tagged or ask your own question.