Nothing
R/filter.C.R
In mskf: Markov Switching Kalman Filter
Defines functions
`filter.C`
`filter.old`
`filter.C` <-
function(theta,npar,ny,ne,nx,nt,nm,mdim,MA,PA,a0,P0,y,x, debug=FALSE)
{
# START WITH MAKING THE ARRAYS THAT ARE NEEDED
auc <- array(0,c((nt + 1),ne,nm)) # a-updated and collapsed to nm posteriors
au <- array(0,c((nt),ne,nm,nm)) # a-updated (not collapsed: S(t-1)=i, S(t)=j)
ap <- array(0,c((nt),ne,nm,nm)) # a-predicted (by definition not collapsed)
Puc <- array(0,c((nt + 1),ne,ne,nm)) # covariance matrix of auc
Pu <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of au
Pp <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of ap
v <- array(0,c(nt,ny,nm,nm)) # one-step-ahead prediction error for S(t-1)=1, S(t)=j
F <- array(0,c(nt,ny,ny,nm,nm)) # covariance matrix of v
pr <- array(0,c((nt + 1),nm)) # updated probabilities
tpr <- array(0,c(nt,nm,nm)) # conditional transition probabilites
utpr <- array(0,c(nt,nm,nm)) # updated conditional transition probabilites
jd <- array(0,c(nt,nm,nm)) # joint densities (product of conditional density and transition probability)
md <- array(0,c(nt,1)) # marginal density (used in optimization)
# MAKE UPDATED PROABILITIES (I.E., VECTOR WITH PROBABILITIES OF
# OF BEING IN CERTAIN STATE AT OCCASION 0)
pr[1,] <- matrix(c(1 / nm),nm,1)
# ENTER INITIAL VALUES OF a0 AND P0 IN auc
auc[1,,] <- a0
Puc[1,,,] <- P0
# PLACE THE ELEMENTS FROM x (= vector with free parameters)
# IN THE APPROPRIATE MODEL MATRICES
for (m in 1:7)
{
for (s in 1:nm)
{
for (r in 1:mdim[m,1])
{
for (c in 1:mdim[m,2])
{
p <- PA[m,r,c,s]
if (p > 0) {
MA[m,r,c,s] <- theta[p]
}
}
}
}
}
# TO OBTAIN THE MATRIX WITH TRANSITION PROBABILITIES,
# WE HAVE TO TAKE INTO ACCOUNT THAT THE PORBS HAVE TO SUM TO 1 PER COLUMN
# HENCE, PER COLUMN THERE ARE nm-1 PROBS TO BE ESTIMATED
# RATHER THAN ESTIMATING THE PROBS WE WILL ESTIMATE A NUMBER
# WHICH CAN RUN FROM -INFINITY TO +INFINITY
tp <- matrix(0, nm, nm)
for (c in 1:(nm - 1))
{
tp[nm,c] <- exp(0)
for (r in 1:(nm - 1))
{
npar <- npar + 1
tp[r,c] <- exp(theta[npar])
}
}
tp[1,nm] <- exp(0)
for (r in 2:nm)
{
npar <- npar + 1
tp[r,nm] <- exp(theta[npar])
}
# THE ZEROS ARE INCLUDED FOR IDENTIFICATION
# TO COME FROM THE ABOVE MATRIX tp TO THE MATRIX WITH ACTUAL
# TRANSITION PROBS p, WE DO THE FOLLOWING
p <- matrix(0,nm,nm)
for (c in 1:nm)
{
stp <- sum(tp[,c])
p[,c] <- tp[,c] / stp
}
# SELECT THE MODEL MATRICES AND NAME THEM
W <- MA[1,1:ny,1:ne,,drop = FALSE]
B <- MA[2,1:ny,1:nx,,drop = FALSE]
R <- MA[3,1:ny,1:ny,,drop = FALSE]
c <- MA[4,1:ne,1,,drop = FALSE]
H <- MA[5,1:ne,1:ne,,drop = FALSE]
G <- MA[6,1:ne,1:ne,,drop = FALSE]
K <- MA[7,1:ne,1:ne,,drop = FALSE]
# START THE KALMAN FILTER & REGIME SWITCHING
L <- 0
L <- Kfilter.timeloop(auc, au, ap,
Puc, Pu, Pp,
v, F, p,
pr, tpr, utpr, jd, md,
W, B, R, c, H, G, K, y, x, debug = debug)
# for (t in 1:nt)
# { # DETERMINE WHETHER THERE ARE MISSING VALUES
# if(all(!y[t,]==-999))
# { for (j in 1:nm)
# { for (i in 1:nm)
# { # KALMAN FILTER STEP
## PREDICTION EQUATIONS
# ap[t,,i,j]<-c[1,,1,j]+H[1,,,j]%*%(auc[t,,i])
# Pp[t,,,i,j]<-H[1,,,j]%*%Puc[t,,,i]%*%t(H[1,,,j])+G[1,,,j]%*%K[1,,,j]%*%t(G[1,,,j])
## ONE-STEP-AHEAD PREDICTION ERROR
# v[t,,i,j]<-y[t,]-W[1,,,j]%*%ap[t,,i,j]-as.matrix(B[1,,,j])%*%x[t,]
# F[t,,,i,j]<-W[1,,,j]%*%Pp[t,,,i,j]%*%t(W[1,,,j])+R[1,,,j]
## UPDATE EQUATIONS
# au[t,,i,j]<-ap[t,,i,j]+Pp[t,,,i,j]%*%t(W[1,,,j])%*%solve(F[t,,,i,j])%*%v[t,,i,j]
# Pu[t,,,i,j]<-Pp[t,,,i,j]+Pp[t,,,i,j]%*%t(W[1,,,j])%*%solve(F[t,,,i,j])%*%W[1,,,j]%*%Pp[t,,,i,j]
## END OF KALMAN FILTER STEP
## BEGIN HAMILTON FILTER STEP
## CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
# tpr[t,j,i]<-p[j,i]*pr[t,i]
## CONDITIONAL DENSITY
# jd[t,j,i]<-tpr[t,j,i]*((2*pi)^(-ny/2)*(det(as.matrix(F[t,,,i,j])))^(-1/2)*exp((-1/2)*t(v[t,,i,j])%*%solve(F[t,,,i,j])%*%v[t,,i,j]))
# }
# }
# jm<-sum(jd[t,,])
# L<-L+log(jm) # LOG LIKELIHOOD FUNCTION
# utpr[t,,]<-jd[t,,]/jm
## END OF THE HAMILTON FILTER STEP
## BEGIN COLLAPSE STEP to reduce the number of posteriors
# for (j in 1:nm)
# { pr[(t+1),j]<-sum(utpr[t,j,])
# as<-matrix(c(0),ne,1)
# for (i in 1:nm)
# { as<-as+au[t,,i,j]*(utpr[t,j,i]/pr[(t+1),j])
# }
# auc[(t+1),,j]<-as
# Ps<-matrix(c(0),ne,ne)
# for (i in 1:nm)
# { Ps<-Ps+(Pu[t,,,i,j]+(auc[(t+1),,j]-au[t,,i,j])%*%(t(auc[(t+1),,j]-au[t,,i,j])))*(utpr[t,j,i]/pr[(t+1),j])
# }
# Puc[(t+1),,,j]<-Ps
# }
# }
# if(any(y[t,]==-999))
# { for (j in 1:nm)
# { for (i in 1:nm)
# { # KALMAN FILTER STEP
## PREDICTION EQUATIONS
# ap[t,,i,j]<-c[1,,1,j]+H[1,,,j]%*%(auc[t,,i])
# Pp[t,,,i,j]<-H[1,,,j]%*%Puc[t,,,i]%*%t(H[1,,,j])+G[1,,,j]%*%K[1,,,j]%*%t(G[1,,,j])
## ONE-STEP-AHEAD PREDICTION ERROR CAN NOT BE CALCULATED
## BUT THE ASSOCIATED COVARIANCE MATRIX CAN
# F[t,,,i,j]<-W[1,,,j]%*%Pp[t,,,i,j]%*%t(W[1,,,j])+R[1,,,j]
## THUS, THE UPDATED STATE IS IDENTICAL TO THE PREDICTED STATE
## AND THE ASSOCIATED COVARIANCE MATRIX IS CALCULATED AS USUAL
# au[t,,i,j]<-ap[t,,i,j]
# Pu[t,,,i,j]<-Pp[t,,,i,j]+Pp[t,,,i,j]%*%t(W[1,,,j])%*%solve(F[t,,,i,j])%*%W[1,,,j]%*%Pp[t,,,i,j]
## END OF KALMAN FILTER STEP
## BEGIN HAMILTON FILTER STEP
## CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
# tpr[t,j,i]<-p[j,i]*pr[t,i]
## CONDITIONAL DENSITY CAN NOT BE CALCULATED
# }
# }
## CURRENT OCCASION DOES NOT CONTRIBUTE TO THE LIKELIHOOD FUNCTION
# utpr[t,,]<-tpr[t,,]
## END OF THE HAMILTON FILTER STEP
## BEGIN COLLAPSE STEP to reduce the number of posteriors
# for (j in 1:nm)
# { pr[(t+1),j]<-sum(utpr[t,j,])
# as<-matrix(c(0),ne,1)
# for (i in 1:nm)
# { as<-as+au[t,,i,j]*(utpr[t,j,i]/pr[(t+1),j])
# }
# auc[(t+1),,j]<-as
# Ps<-matrix(c(0),ne,ne)
# for (i in 1:nm)
# { Ps<-Ps+(Pu[t,,,i,j]+(auc[(t+1),,j]-au[t,,i,j])%*%(t(auc[(t+1),,j]-au[t,,i,j])))*(utpr[t,j,i]/pr[(t+1),j])
# }
# Puc[(t+1),,,j]<-Ps
# }
# }
# } # END OF TIME LOOP
Lp<-if(is.nan(L)) 1e9 else min(-L,1e9)
attr(Lp,"jd") <- jd
attr(Lp,"all") <- list(jd=jd, Puc=Puc, Pu=Pu, Pp=Pp, v=v, F=F, pr=pr, tpr=tpr, utpr=utpr, au=au, ap=ap, auc=auc, md=md)
Lp
}
`filter.old` <-
function(theta,npar,ny,ne,nx,nt,nm,mdim,MA,PA,a0,P0,y,x, debug=FALSE)
{
Debug = debug
# START WITH MAKING THE ARRAYS THAT ARE NEEDED
auc <- array(0,c((nt + 1),ne,nm)) # a-updated and collapsed to nm posteriors
au <- array(0,c((nt),ne,nm,nm)) # a-updated (not collapsed: S(t-1)=i, S(t)=j)
ap <- array(0,c((nt),ne,nm,nm)) # a-predicted (by definition not collapsed)
Puc <- array(0,c((nt + 1),ne,ne,nm)) # covariance matrix of auc
Pu <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of au
Pp <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of ap
v <- array(0,c(nt,ny,nm,nm)) # one-step-ahead prediction error for S(t-1)=1, S(t)=j
F <- array(0,c(nt,ny,ny,nm,nm)) # covariance matrix of v
pr <- array(0,c((nt + 1),nm)) # updated probabilities
tpr <- array(0,c(nt,nm,nm)) # conditional transition probabilites
utpr <- array(0,c(nt,nm,nm)) # updated conditional transition probabilites
jd <- array(0,c(nt,nm,nm)) # joint densities (product of conditional density and transition probability)
md <- array(0,c(nt,1)) # marginal density (used in optimization)
# MAKE UPDATED PROABILITIES (I.E., VECTOR WITH PROBABILITIES OF
# OF BEING IN CERTAIN STATE AT OCCASION 0)
pr[1,] <- matrix(c(1 / nm),nm,1)
# ENTER INITIAL VALUES OF a0 AND P0 IN auc
auc[1,,] <- a0
Puc[1,,,] <- P0
# PLACE THE ELEMENTS FROM x (= vector with free parameters)
# IN THE APPROPRIATE MODEL MATRICES
for (m in 1:7)
{
for (s in 1:nm)
{
for (r in 1:mdim[m,1])
{
for (c in 1:mdim[m,2])
{
p <- PA[m,r,c,s]
if (p > 0) {
MA[m,r,c,s] <- theta[p]
}
}
}
}
}
# TO OBTAIN THE MATRIX WITH TRANSITION PROBABILITIES,
# WE HAVE TO TAKE INTO ACCOUNT THAT THE PORBS HAVE TO SUM TO 1 PER COLUMN
# HENCE, PER COLUMN THERE ARE nm-1 PROBS TO BE ESTIMATED
# RATHER THAN ESTIMATING THE PROBS WE WILL ESTIMATE A NUMBER
# WHICH CAN RUN FROM -INFINITY TO +INFINITY
tp <- matrix(0, nm, nm)
for (c in 1:(nm - 1))
{
tp[nm,c] <- exp(0)
for (r in 1:(nm - 1))
{
npar <- npar + 1
tp[r,c] <- exp(theta[npar])
}
}
tp[1,nm] <- exp(0)
for (r in 2:nm)
{
npar <- npar + 1
tp[r,nm] <- exp(theta[npar])
}
# THE ZEROS ARE INCLUDED FOR IDENTIFICATION
# TO COME FROM THE ABOVE MATRIX tp TO THE MATRIX WITH ACTUAL
# TRANSITION PROBS p, WE DO THE FOLLOWING
p <- matrix(0,nm,nm)
for (c in 1:nm)
{
stp <- sum(tp[,c])
p[,c] <- tp[,c] / stp
}
# SELECT THE MODEL MATRICES AND NAME THEM
W <- MA[1,1:ny,1:ne,,drop = FALSE]
B <- MA[2,1:ny,1:nx,,drop = FALSE]
R <- MA[3,1:ny,1:ny,,drop = FALSE]
c <- MA[4,1:ne,1,,drop = FALSE]
H <- MA[5,1:ne,1:ne,,drop = FALSE]
G <- MA[6,1:ne,1:ne,,drop = FALSE]
K <- MA[7,1:ne,1:ne,,drop = FALSE]
# START THE KALMAN FILTER & REGIME SWITCHING
L <- 0
# L <- Kfilter.timeloop(auc, au, ap,
# Puc, Pu, Pp,
# v, F, p,
# pr, tpr, utpr, jd, md,
# W, B, R, c, H, G, K, y, x, debug = debug)
### ----------------- debug code ------------------- ###
if (debug) {
cat("x[0,1] = ", x[0,1], "\n")
cat(paste(c("nt","ne","nm","ny","nx"),c(nt,ne,nm,ny,nx),sep=" = ", collapse = "\n"))
cat("\np =\n"); print(p)
for (j in 1:nm) {
cat(sprintf("\n\nRegime %d\n========\n\n\n\ta[t] = c + H a[t-1] + G u[t], u[t] ~ N(0, K)\n\ty[t] = W a[t] + B x[t] + e[t], e[t] ~ N(0, R)\n\n", j-1));
cat("\nH = \n"); print(as.matrix(H[1,,,j]));
cat("\nK =\n"); print(as.matrix(K[1,,,j]));
cat("\nG =\n"); print(as.matrix(G[1,,,j]));
cat("\nc =\n"); print(as.matrix(c[1,,,j]));
cat("\nB =\n"); print(as.matrix(B[1,,,j]));
cat("\nR =\n"); print(as.matrix(R[1,,,j]));
cat("\nW =\n"); print(as.matrix(W[1,,,j]));
}
cat("-----------------------------------------------------------------------------------------------------------\n");
Debug = debug;
}
### ----------------- end debug ------------------- ###
for (t in 1:nt)
{
### ----------------- debug code ------------------- ###
if (t >= debug) debug = 0;
if (debug) {
cat(sprintf("\nt = %d\n", t-1));
}
### ----------------- end debug ------------------- ###
# DETERMINE WHETHER THERE ARE MISSING VALUES
if (all(!is.na(y[t,])))
{
for (j in 1:nm)
{
for (i in 1:nm)
{
# KALMAN FILTER STEP
# PREDICTION EQUATIONS
### ----------------- debug code ------------------- ###
if(debug && i==1 && j==1) {
cat("\np = \n"); print(p);
}
### ----------------- end debug ------------------- ###
ap[t,,i,j] <- c[1,,1,j] + H[1,,,j] %*% (auc[t,,i])
### ----------------- debug code ------------------- ###
if(debug) {
cat("\nap = c + H %*% auc\n");
cat("\nc = \n"); print(as.matrix(c[1,,,j]));
cat("\nH = \n"); print(as.matrix(H[1,,,j]));
cat("\nauc = \n"); print(as.matrix(auc[t,,j]));
cat("\nap = \n"); print(as.matrix(ap[t,,i,j]));
}
### ----------------- end debug ------------------- ###
Pp[t,,,i,j] <-
H[1,,,j] %*% Puc[t,,,i] %*% t(H[1,,,j]) + G[1,,,j] %*% K[1,,,j] %*% t(G[1,,,j])
### ----------------- debug code ------------------- ###
if(debug) cat(sprintf("\n S(i=%d) --> S(j=%d) at t = %d\t\t***********************************************\n\nPp = H Puc H' + G K G'", i-1, j-1, t-1)); # states are labeled 0, 1, 2, etc
if(debug){
cat("\n");
cat("\nPuc = \n"); print(as.matrix(Puc[t,,,j]));
cat("\nG = \n"); print(as.matrix(G[1,,,j]));
cat("\nK = \n"); print(as.matrix(K[1,,,j]));
cat("\nPp = \n"); print(as.matrix(Pp[t,,,i,j]));
}
### ----------------- end debug ------------------- ###
# ONE-STEP-AHEAD PREDICTION ERROR
v[t,,i,j] <- y[t,] - W[1,,,j] %*% ap[t,,i,j] - as.matrix(B[1,,,j]) %*% x[t,]
F[t,,,i,j] <- W[1,,,j] %*% Pp[t,,,i,j] %*% t(W[1,,,j]) + R[1,,,j]
### ----------------- debug code ------------------- ###
if(debug){
cat("\nv1 = y - W ap\n");
cat("\ny = \n"); print(as.matrix(y[t,]))
cat("\nW = \n"); print(as.matrix(W[1,,,j]));
cat("\nap = \n"); print(as.matrix(ap[t,,i,j]));
cat("\nv1 = \n"); print(matrix(y[t,] - W[1,,,j] %*% ap[t,,i,j]));
cat("\nv = v1 - B x\n");
cat("\nB = \n"); print(as.matrix(B[1,,,j]));
cat("\nx = \n"); print(as.matrix(x[t,]))
cat("\nv = \n"); print(as.matrix(v[t,,i,j]));
cat("\nF = W Pp W' + R\n");
cat("\nPp = \n"); print(as.matrix(Pp[t,,,i,j]));
cat("\nPpW = \n"); print(as.matrix(Pp[t,,,i,j] %*% t(W[1,,,j])));
cat("\nR = \n"); print(as.matrix(R[1,,,j]));
cat("\nF = \n"); print(as.matrix(F[t,,,i,j]));
}
### ----------------- end debug ------------------- ###
# UPDATE EQUATIONS
au[t,,i,j] <-
ap[t,,i,j] + Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j]
Pu[t,,,i,j] <-
Pp[t,,,i,j] + Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j]) %*% W[1,,,j] %*% Pp[t,,,i,j]
### ----------------- debug code ------------------- ###
if(debug){
cat("\n");
cat("\nFinv =\n"); print(as.matrix(solve(F[t,,,i,j])))
cat(sprintf("\n|F| = %12.6g\n", det(F[t,,,i,j])));
cat("\nau = ap + Pp W Finv v\n");
cat("\nPpWFinv = \n"); print(as.matrix(Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j])));
cat("\nau = \n"); print(as.matrix(au[t, , i, j]));
cat("Pu = Pp - Pp W Finv W' Pp\n");
cat("\nPu = \n"); print(as.matrix(Pu[t,,,i,j]));
}
### ----------------- end debug ------------------- ###
# END OF KALMAN FILTER STEP
# BEGIN HAMILTON FILTER STEP
# CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
tpr[t,j,i] <- p[j,i] * pr[t,i]
# CONDITIONAL DENSITY
jd[t,j,i] <-
tpr[t,j,i] * ((2 * pi) ^ (-ny / 2) * det(as.matrix(F[t,,,i,j]))^(-1/2) * exp((-1/2) * t(v[t,,i,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j]))
### ----------------- debug code ------------------- ***/
if (Debug) {
cat(sprintf("\n v = %12.6g\nv' Finv v = %12.6g\n", v[t,1,i,j], t(v[t,,i,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j]));
cat(sprintf("\njd[%d,%d] = tpr[%d,%d] * (2 * pi)^(ny/2) * exp(-0.5 * v' Finv v) / sqrt(det(F))", j, i, j, i));
cat(sprintf("\n\t p = %12.6g\n\t pr = %12.6g\n\t tpr = %12.6g\n\t(2*pi)^(-ny/2) = %12.6g\n\t vtFinvv[0] = %12.6g\n\t detF = %12.6g\n\t jd[%d,%d] = %12.6g\n", p[j, i], pr[t, i], tpr[t,j,i], (2*pi)^(-ny/2.0), t(v[t,,i,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j], det(as.matrix(F[t,,,i,j])), j, i, jd[t, j, i]));
}
### ----------------- end debug ------------------- ***/
}
}
jm <- sum(jd[t,,])
L <- L + log(jm) # LOG LIKELIHOOD FUNCTION
utpr[t,,] <- jd[t,,] / jm
### ----------------- debug code ------------------- ***/
if (debug) {
cat(sprintf("\n\n +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\t missing = %12s\n\tlog-Likelihood = %12.6g\n\tjm = %12.6g\n\n", if (any(is.na(y[t,]))) "TRUE" else "FALSE", L, jm));
}
### ----------------- end debug ------------------- ***/
# END OF THE HAMILTON FILTER STEP
# BEGIN COLLAPSE STEP to reduce the number of posteriors
for (j in 1:nm)
{
pr[(t + 1),j] <- sum(utpr[t,j,])
as <- matrix(c(0),ne,1)
for (i in 1:nm)
{
as <- as + au[t,,i,j] * (utpr[t,j,i] / pr[(t + 1),j])
}
auc[(t + 1),,j] <- as
Ps <- matrix(c(0),ne,ne)
for (i in 1:nm)
{
Ps <-
Ps + (Pu[t,,,i,j] + (auc[(t + 1),,j] - au[t,,i,j]) %*% (t(auc[(t + 1),,j] -
au[t,,i,j]))) * (utpr[t,j,i] / pr[(t + 1),j])
}
Puc[(t + 1),,,j] <- Ps
}
}
###
### If there are missing values we have a separate time loop... (for no reason)
###
if (any(is.na(y[t,])))
{
### ----------------- debug code ------------------- ###
if(debug)
cat(sprintf("missing value(s) at t = %d\n", t));
### ----------------- end debug ------------------- ###
for (j in 1:nm)
{
for (i in 1:nm)
{
# KALMAN FILTER STEP
# PREDICTION EQUATIONS
ap[t,,i,j] <- c[1,,1,j] + H[1,,,j] %*% (auc[t,,i])
Pp[t,,,i,j] <-
H[1,,,j] %*% Puc[t,,,i] %*% t(H[1,,,j]) + G[1,,,j] %*% K[1,,,j] %*% t(G[1,,,j])
# ONE-STEP-AHEAD PREDICTION ERROR CAN NOT BE CALCULATED
# BUT THE ASSOCIATED COVARIANCE MATRIX CAN
F[t,,,i,j] <- W[1,,,j] %*% Pp[t,,,i,j] %*% t(W[1,,,j]) + R[1,,,j]
# THUS, THE UPDATED STATE IS IDENTICAL TO THE PREDICTED STATE
# AND THE ASSOCIATED COVARIANCE MATRIX IS CALCULATED AS USUAL
au[t,,i,j] <- ap[t,,i,j]
Pu[t,,,i,j] <-
Pp[t,,,i,j] + Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j]) %*% W[1,,,j] %*% Pp[t,,,i,j]
# END OF KALMAN FILTER STEP
# BEGIN HAMILTON FILTER STEP
# CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
tpr[t,j,i] <- p[j,i] * pr[t,i]
# CONDITIONAL DENSITY CAN NOT BE CALCULATED
}
}
# CURRENT OCCASION DOES NOT CONTRIBUTE TO THE LIKELIHOOD FUNCTION
utpr[t,,] <- tpr[t,,]
# END OF THE HAMILTON FILTER STEP
# BEGIN COLLAPSE STEP to reduce the number of posteriors
for (j in 1:nm)
{
pr[(t + 1),j] <- sum(utpr[t,j,])
as <- matrix(c(0),ne,1)
for (i in 1:nm)
{
as <- as + au[t,,i,j] * (utpr[t,j,i] / pr[(t + 1),j])
}
auc[(t + 1),,j] <- as
Ps <- matrix(c(0),ne,ne)
for (i in 1:nm)
{
Ps <-
Ps + (Pu[t,,,i,j] + (auc[(t + 1),,j] - au[t,,i,j]) %*% (t(auc[(t + 1),,j] -
au[t,,i,j]))) * (utpr[t,j,i] / pr[(t + 1),j])
}
Puc[(t + 1),,,j] <- Ps
}
}
} # END OF TIME LOOP
Lp<-if(is.nan(L)) 1e9 else min(-L,1e9)
attr(Lp,"jd") <- jd
attr(Lp,"all") <- list(jd=jd, Puc=Puc, Pu=Pu, Pp=Pp, v=v, F=F, pr=pr, tpr=tpr, utpr=utpr, au=au, ap=ap, auc=auc, md=md)
Lp
}
Try the mskf package in your browser
Any scripts or data that you put into this service are public.
mskf documentation built on May 2, 2019, 6:47 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions `filter.C` `filter.old`
`filter.C` <-
function(theta,npar,ny,ne,nx,nt,nm,mdim,MA,PA,a0,P0,y,x, debug=FALSE)
{
# START WITH MAKING THE ARRAYS THAT ARE NEEDED
auc <- array(0,c((nt + 1),ne,nm)) # a-updated and collapsed to nm posteriors
au <- array(0,c((nt),ne,nm,nm)) # a-updated (not collapsed: S(t-1)=i, S(t)=j)
ap <- array(0,c((nt),ne,nm,nm)) # a-predicted (by definition not collapsed)
Puc <- array(0,c((nt + 1),ne,ne,nm)) # covariance matrix of auc
Pu <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of au
Pp <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of ap
v <- array(0,c(nt,ny,nm,nm)) # one-step-ahead prediction error for S(t-1)=1, S(t)=j
F <- array(0,c(nt,ny,ny,nm,nm)) # covariance matrix of v
pr <- array(0,c((nt + 1),nm)) # updated probabilities
tpr <- array(0,c(nt,nm,nm)) # conditional transition probabilites
utpr <- array(0,c(nt,nm,nm)) # updated conditional transition probabilites
jd <- array(0,c(nt,nm,nm)) # joint densities (product of conditional density and transition probability)
md <- array(0,c(nt,1)) # marginal density (used in optimization)
# MAKE UPDATED PROABILITIES (I.E., VECTOR WITH PROBABILITIES OF
# OF BEING IN CERTAIN STATE AT OCCASION 0)
pr[1,] <- matrix(c(1 / nm),nm,1)
# ENTER INITIAL VALUES OF a0 AND P0 IN auc
auc[1,,] <- a0
Puc[1,,,] <- P0
# PLACE THE ELEMENTS FROM x (= vector with free parameters)
# IN THE APPROPRIATE MODEL MATRICES
for (m in 1:7)
{
for (s in 1:nm)
{
for (r in 1:mdim[m,1])
{
for (c in 1:mdim[m,2])
{
p <- PA[m,r,c,s]
if (p > 0) {
MA[m,r,c,s] <- theta[p]
}
}
}
}
}
# TO OBTAIN THE MATRIX WITH TRANSITION PROBABILITIES,
# WE HAVE TO TAKE INTO ACCOUNT THAT THE PORBS HAVE TO SUM TO 1 PER COLUMN
# HENCE, PER COLUMN THERE ARE nm-1 PROBS TO BE ESTIMATED
# RATHER THAN ESTIMATING THE PROBS WE WILL ESTIMATE A NUMBER
# WHICH CAN RUN FROM -INFINITY TO +INFINITY
tp <- matrix(0, nm, nm)
for (c in 1:(nm - 1))
{
tp[nm,c] <- exp(0)
for (r in 1:(nm - 1))
{
npar <- npar + 1
tp[r,c] <- exp(theta[npar])
}
}
tp[1,nm] <- exp(0)
for (r in 2:nm)
{
npar <- npar + 1
tp[r,nm] <- exp(theta[npar])
}
# THE ZEROS ARE INCLUDED FOR IDENTIFICATION
# TO COME FROM THE ABOVE MATRIX tp TO THE MATRIX WITH ACTUAL
# TRANSITION PROBS p, WE DO THE FOLLOWING
p <- matrix(0,nm,nm)
for (c in 1:nm)
{
stp <- sum(tp[,c])
p[,c] <- tp[,c] / stp
}
# SELECT THE MODEL MATRICES AND NAME THEM
W <- MA[1,1:ny,1:ne,,drop = FALSE]
B <- MA[2,1:ny,1:nx,,drop = FALSE]
R <- MA[3,1:ny,1:ny,,drop = FALSE]
c <- MA[4,1:ne,1,,drop = FALSE]
H <- MA[5,1:ne,1:ne,,drop = FALSE]
G <- MA[6,1:ne,1:ne,,drop = FALSE]
K <- MA[7,1:ne,1:ne,,drop = FALSE]
# START THE KALMAN FILTER & REGIME SWITCHING
L <- 0
L <- Kfilter.timeloop(auc, au, ap,
Puc, Pu, Pp,
v, F, p,
pr, tpr, utpr, jd, md,
W, B, R, c, H, G, K, y, x, debug = debug)
# for (t in 1:nt)
# { # DETERMINE WHETHER THERE ARE MISSING VALUES
# if(all(!y[t,]==-999))
# { for (j in 1:nm)
# { for (i in 1:nm)
# { # KALMAN FILTER STEP
## PREDICTION EQUATIONS
# ap[t,,i,j]<-c[1,,1,j]+H[1,,,j]%*%(auc[t,,i])
# Pp[t,,,i,j]<-H[1,,,j]%*%Puc[t,,,i]%*%t(H[1,,,j])+G[1,,,j]%*%K[1,,,j]%*%t(G[1,,,j])
## ONE-STEP-AHEAD PREDICTION ERROR
# v[t,,i,j]<-y[t,]-W[1,,,j]%*%ap[t,,i,j]-as.matrix(B[1,,,j])%*%x[t,]
# F[t,,,i,j]<-W[1,,,j]%*%Pp[t,,,i,j]%*%t(W[1,,,j])+R[1,,,j]
## UPDATE EQUATIONS
# au[t,,i,j]<-ap[t,,i,j]+Pp[t,,,i,j]%*%t(W[1,,,j])%*%solve(F[t,,,i,j])%*%v[t,,i,j]
# Pu[t,,,i,j]<-Pp[t,,,i,j]+Pp[t,,,i,j]%*%t(W[1,,,j])%*%solve(F[t,,,i,j])%*%W[1,,,j]%*%Pp[t,,,i,j]
## END OF KALMAN FILTER STEP
## BEGIN HAMILTON FILTER STEP
## CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
# tpr[t,j,i]<-p[j,i]*pr[t,i]
## CONDITIONAL DENSITY
# jd[t,j,i]<-tpr[t,j,i]*((2*pi)^(-ny/2)*(det(as.matrix(F[t,,,i,j])))^(-1/2)*exp((-1/2)*t(v[t,,i,j])%*%solve(F[t,,,i,j])%*%v[t,,i,j]))
# }
# }
# jm<-sum(jd[t,,])
# L<-L+log(jm) # LOG LIKELIHOOD FUNCTION
# utpr[t,,]<-jd[t,,]/jm
## END OF THE HAMILTON FILTER STEP
## BEGIN COLLAPSE STEP to reduce the number of posteriors
# for (j in 1:nm)
# { pr[(t+1),j]<-sum(utpr[t,j,])
# as<-matrix(c(0),ne,1)
# for (i in 1:nm)
# { as<-as+au[t,,i,j]*(utpr[t,j,i]/pr[(t+1),j])
# }
# auc[(t+1),,j]<-as
# Ps<-matrix(c(0),ne,ne)
# for (i in 1:nm)
# { Ps<-Ps+(Pu[t,,,i,j]+(auc[(t+1),,j]-au[t,,i,j])%*%(t(auc[(t+1),,j]-au[t,,i,j])))*(utpr[t,j,i]/pr[(t+1),j])
# }
# Puc[(t+1),,,j]<-Ps
# }
# }
# if(any(y[t,]==-999))
# { for (j in 1:nm)
# { for (i in 1:nm)
# { # KALMAN FILTER STEP
## PREDICTION EQUATIONS
# ap[t,,i,j]<-c[1,,1,j]+H[1,,,j]%*%(auc[t,,i])
# Pp[t,,,i,j]<-H[1,,,j]%*%Puc[t,,,i]%*%t(H[1,,,j])+G[1,,,j]%*%K[1,,,j]%*%t(G[1,,,j])
## ONE-STEP-AHEAD PREDICTION ERROR CAN NOT BE CALCULATED
## BUT THE ASSOCIATED COVARIANCE MATRIX CAN
# F[t,,,i,j]<-W[1,,,j]%*%Pp[t,,,i,j]%*%t(W[1,,,j])+R[1,,,j]
## THUS, THE UPDATED STATE IS IDENTICAL TO THE PREDICTED STATE
## AND THE ASSOCIATED COVARIANCE MATRIX IS CALCULATED AS USUAL
# au[t,,i,j]<-ap[t,,i,j]
# Pu[t,,,i,j]<-Pp[t,,,i,j]+Pp[t,,,i,j]%*%t(W[1,,,j])%*%solve(F[t,,,i,j])%*%W[1,,,j]%*%Pp[t,,,i,j]
## END OF KALMAN FILTER STEP
## BEGIN HAMILTON FILTER STEP
## CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
# tpr[t,j,i]<-p[j,i]*pr[t,i]
## CONDITIONAL DENSITY CAN NOT BE CALCULATED
# }
# }
## CURRENT OCCASION DOES NOT CONTRIBUTE TO THE LIKELIHOOD FUNCTION
# utpr[t,,]<-tpr[t,,]
## END OF THE HAMILTON FILTER STEP
## BEGIN COLLAPSE STEP to reduce the number of posteriors
# for (j in 1:nm)
# { pr[(t+1),j]<-sum(utpr[t,j,])
# as<-matrix(c(0),ne,1)
# for (i in 1:nm)
# { as<-as+au[t,,i,j]*(utpr[t,j,i]/pr[(t+1),j])
# }
# auc[(t+1),,j]<-as
# Ps<-matrix(c(0),ne,ne)
# for (i in 1:nm)
# { Ps<-Ps+(Pu[t,,,i,j]+(auc[(t+1),,j]-au[t,,i,j])%*%(t(auc[(t+1),,j]-au[t,,i,j])))*(utpr[t,j,i]/pr[(t+1),j])
# }
# Puc[(t+1),,,j]<-Ps
# }
# }
# } # END OF TIME LOOP
Lp<-if(is.nan(L)) 1e9 else min(-L,1e9)
attr(Lp,"jd") <- jd
attr(Lp,"all") <- list(jd=jd, Puc=Puc, Pu=Pu, Pp=Pp, v=v, F=F, pr=pr, tpr=tpr, utpr=utpr, au=au, ap=ap, auc=auc, md=md)
Lp
}
`filter.old` <-
function(theta,npar,ny,ne,nx,nt,nm,mdim,MA,PA,a0,P0,y,x, debug=FALSE)
{
Debug = debug
# START WITH MAKING THE ARRAYS THAT ARE NEEDED
auc <- array(0,c((nt + 1),ne,nm)) # a-updated and collapsed to nm posteriors
au <- array(0,c((nt),ne,nm,nm)) # a-updated (not collapsed: S(t-1)=i, S(t)=j)
ap <- array(0,c((nt),ne,nm,nm)) # a-predicted (by definition not collapsed)
Puc <- array(0,c((nt + 1),ne,ne,nm)) # covariance matrix of auc
Pu <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of au
Pp <- array(0,c(nt,ne,ne,nm,nm)) # covariance matrix of ap
v <- array(0,c(nt,ny,nm,nm)) # one-step-ahead prediction error for S(t-1)=1, S(t)=j
F <- array(0,c(nt,ny,ny,nm,nm)) # covariance matrix of v
pr <- array(0,c((nt + 1),nm)) # updated probabilities
tpr <- array(0,c(nt,nm,nm)) # conditional transition probabilites
utpr <- array(0,c(nt,nm,nm)) # updated conditional transition probabilites
jd <- array(0,c(nt,nm,nm)) # joint densities (product of conditional density and transition probability)
md <- array(0,c(nt,1)) # marginal density (used in optimization)
# MAKE UPDATED PROABILITIES (I.E., VECTOR WITH PROBABILITIES OF
# OF BEING IN CERTAIN STATE AT OCCASION 0)
pr[1,] <- matrix(c(1 / nm),nm,1)
# ENTER INITIAL VALUES OF a0 AND P0 IN auc
auc[1,,] <- a0
Puc[1,,,] <- P0
# PLACE THE ELEMENTS FROM x (= vector with free parameters)
# IN THE APPROPRIATE MODEL MATRICES
for (m in 1:7)
{
for (s in 1:nm)
{
for (r in 1:mdim[m,1])
{
for (c in 1:mdim[m,2])
{
p <- PA[m,r,c,s]
if (p > 0) {
MA[m,r,c,s] <- theta[p]
}
}
}
}
}
# TO OBTAIN THE MATRIX WITH TRANSITION PROBABILITIES,
# WE HAVE TO TAKE INTO ACCOUNT THAT THE PORBS HAVE TO SUM TO 1 PER COLUMN
# HENCE, PER COLUMN THERE ARE nm-1 PROBS TO BE ESTIMATED
# RATHER THAN ESTIMATING THE PROBS WE WILL ESTIMATE A NUMBER
# WHICH CAN RUN FROM -INFINITY TO +INFINITY
tp <- matrix(0, nm, nm)
for (c in 1:(nm - 1))
{
tp[nm,c] <- exp(0)
for (r in 1:(nm - 1))
{
npar <- npar + 1
tp[r,c] <- exp(theta[npar])
}
}
tp[1,nm] <- exp(0)
for (r in 2:nm)
{
npar <- npar + 1
tp[r,nm] <- exp(theta[npar])
}
# THE ZEROS ARE INCLUDED FOR IDENTIFICATION
# TO COME FROM THE ABOVE MATRIX tp TO THE MATRIX WITH ACTUAL
# TRANSITION PROBS p, WE DO THE FOLLOWING
p <- matrix(0,nm,nm)
for (c in 1:nm)
{
stp <- sum(tp[,c])
p[,c] <- tp[,c] / stp
}
# SELECT THE MODEL MATRICES AND NAME THEM
W <- MA[1,1:ny,1:ne,,drop = FALSE]
B <- MA[2,1:ny,1:nx,,drop = FALSE]
R <- MA[3,1:ny,1:ny,,drop = FALSE]
c <- MA[4,1:ne,1,,drop = FALSE]
H <- MA[5,1:ne,1:ne,,drop = FALSE]
G <- MA[6,1:ne,1:ne,,drop = FALSE]
K <- MA[7,1:ne,1:ne,,drop = FALSE]
# START THE KALMAN FILTER & REGIME SWITCHING
L <- 0
# L <- Kfilter.timeloop(auc, au, ap,
# Puc, Pu, Pp,
# v, F, p,
# pr, tpr, utpr, jd, md,
# W, B, R, c, H, G, K, y, x, debug = debug)
### ----------------- debug code ------------------- ###
if (debug) {
cat("x[0,1] = ", x[0,1], "\n")
cat(paste(c("nt","ne","nm","ny","nx"),c(nt,ne,nm,ny,nx),sep=" = ", collapse = "\n"))
cat("\np =\n"); print(p)
for (j in 1:nm) {
cat(sprintf("\n\nRegime %d\n========\n\n\n\ta[t] = c + H a[t-1] + G u[t], u[t] ~ N(0, K)\n\ty[t] = W a[t] + B x[t] + e[t], e[t] ~ N(0, R)\n\n", j-1));
cat("\nH = \n"); print(as.matrix(H[1,,,j]));
cat("\nK =\n"); print(as.matrix(K[1,,,j]));
cat("\nG =\n"); print(as.matrix(G[1,,,j]));
cat("\nc =\n"); print(as.matrix(c[1,,,j]));
cat("\nB =\n"); print(as.matrix(B[1,,,j]));
cat("\nR =\n"); print(as.matrix(R[1,,,j]));
cat("\nW =\n"); print(as.matrix(W[1,,,j]));
}
cat("-----------------------------------------------------------------------------------------------------------\n");
Debug = debug;
}
### ----------------- end debug ------------------- ###
for (t in 1:nt)
{
### ----------------- debug code ------------------- ###
if (t >= debug) debug = 0;
if (debug) {
cat(sprintf("\nt = %d\n", t-1));
}
### ----------------- end debug ------------------- ###
# DETERMINE WHETHER THERE ARE MISSING VALUES
if (all(!is.na(y[t,])))
{
for (j in 1:nm)
{
for (i in 1:nm)
{
# KALMAN FILTER STEP
# PREDICTION EQUATIONS
### ----------------- debug code ------------------- ###
if(debug && i==1 && j==1) {
cat("\np = \n"); print(p);
}
### ----------------- end debug ------------------- ###
ap[t,,i,j] <- c[1,,1,j] + H[1,,,j] %*% (auc[t,,i])
### ----------------- debug code ------------------- ###
if(debug) {
cat("\nap = c + H %*% auc\n");
cat("\nc = \n"); print(as.matrix(c[1,,,j]));
cat("\nH = \n"); print(as.matrix(H[1,,,j]));
cat("\nauc = \n"); print(as.matrix(auc[t,,j]));
cat("\nap = \n"); print(as.matrix(ap[t,,i,j]));
}
### ----------------- end debug ------------------- ###
Pp[t,,,i,j] <-
H[1,,,j] %*% Puc[t,,,i] %*% t(H[1,,,j]) + G[1,,,j] %*% K[1,,,j] %*% t(G[1,,,j])
### ----------------- debug code ------------------- ###
if(debug) cat(sprintf("\n S(i=%d) --> S(j=%d) at t = %d\t\t***********************************************\n\nPp = H Puc H' + G K G'", i-1, j-1, t-1)); # states are labeled 0, 1, 2, etc
if(debug){
cat("\n");
cat("\nPuc = \n"); print(as.matrix(Puc[t,,,j]));
cat("\nG = \n"); print(as.matrix(G[1,,,j]));
cat("\nK = \n"); print(as.matrix(K[1,,,j]));
cat("\nPp = \n"); print(as.matrix(Pp[t,,,i,j]));
}
### ----------------- end debug ------------------- ###
# ONE-STEP-AHEAD PREDICTION ERROR
v[t,,i,j] <- y[t,] - W[1,,,j] %*% ap[t,,i,j] - as.matrix(B[1,,,j]) %*% x[t,]
F[t,,,i,j] <- W[1,,,j] %*% Pp[t,,,i,j] %*% t(W[1,,,j]) + R[1,,,j]
### ----------------- debug code ------------------- ###
if(debug){
cat("\nv1 = y - W ap\n");
cat("\ny = \n"); print(as.matrix(y[t,]))
cat("\nW = \n"); print(as.matrix(W[1,,,j]));
cat("\nap = \n"); print(as.matrix(ap[t,,i,j]));
cat("\nv1 = \n"); print(matrix(y[t,] - W[1,,,j] %*% ap[t,,i,j]));
cat("\nv = v1 - B x\n");
cat("\nB = \n"); print(as.matrix(B[1,,,j]));
cat("\nx = \n"); print(as.matrix(x[t,]))
cat("\nv = \n"); print(as.matrix(v[t,,i,j]));
cat("\nF = W Pp W' + R\n");
cat("\nPp = \n"); print(as.matrix(Pp[t,,,i,j]));
cat("\nPpW = \n"); print(as.matrix(Pp[t,,,i,j] %*% t(W[1,,,j])));
cat("\nR = \n"); print(as.matrix(R[1,,,j]));
cat("\nF = \n"); print(as.matrix(F[t,,,i,j]));
}
### ----------------- end debug ------------------- ###
# UPDATE EQUATIONS
au[t,,i,j] <-
ap[t,,i,j] + Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j]
Pu[t,,,i,j] <-
Pp[t,,,i,j] + Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j]) %*% W[1,,,j] %*% Pp[t,,,i,j]
### ----------------- debug code ------------------- ###
if(debug){
cat("\n");
cat("\nFinv =\n"); print(as.matrix(solve(F[t,,,i,j])))
cat(sprintf("\n|F| = %12.6g\n", det(F[t,,,i,j])));
cat("\nau = ap + Pp W Finv v\n");
cat("\nPpWFinv = \n"); print(as.matrix(Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j])));
cat("\nau = \n"); print(as.matrix(au[t, , i, j]));
cat("Pu = Pp - Pp W Finv W' Pp\n");
cat("\nPu = \n"); print(as.matrix(Pu[t,,,i,j]));
}
### ----------------- end debug ------------------- ###
# END OF KALMAN FILTER STEP
# BEGIN HAMILTON FILTER STEP
# CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
tpr[t,j,i] <- p[j,i] * pr[t,i]
# CONDITIONAL DENSITY
jd[t,j,i] <-
tpr[t,j,i] * ((2 * pi) ^ (-ny / 2) * det(as.matrix(F[t,,,i,j]))^(-1/2) * exp((-1/2) * t(v[t,,i,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j]))
### ----------------- debug code ------------------- ***/
if (Debug) {
cat(sprintf("\n v = %12.6g\nv' Finv v = %12.6g\n", v[t,1,i,j], t(v[t,,i,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j]));
cat(sprintf("\njd[%d,%d] = tpr[%d,%d] * (2 * pi)^(ny/2) * exp(-0.5 * v' Finv v) / sqrt(det(F))", j, i, j, i));
cat(sprintf("\n\t p = %12.6g\n\t pr = %12.6g\n\t tpr = %12.6g\n\t(2*pi)^(-ny/2) = %12.6g\n\t vtFinvv[0] = %12.6g\n\t detF = %12.6g\n\t jd[%d,%d] = %12.6g\n", p[j, i], pr[t, i], tpr[t,j,i], (2*pi)^(-ny/2.0), t(v[t,,i,j]) %*% solve(F[t,,,i,j]) %*% v[t,,i,j], det(as.matrix(F[t,,,i,j])), j, i, jd[t, j, i]));
}
### ----------------- end debug ------------------- ***/
}
}
jm <- sum(jd[t,,])
L <- L + log(jm) # LOG LIKELIHOOD FUNCTION
utpr[t,,] <- jd[t,,] / jm
### ----------------- debug code ------------------- ***/
if (debug) {
cat(sprintf("\n\n +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\t missing = %12s\n\tlog-Likelihood = %12.6g\n\tjm = %12.6g\n\n", if (any(is.na(y[t,]))) "TRUE" else "FALSE", L, jm));
}
### ----------------- end debug ------------------- ***/
# END OF THE HAMILTON FILTER STEP
# BEGIN COLLAPSE STEP to reduce the number of posteriors
for (j in 1:nm)
{
pr[(t + 1),j] <- sum(utpr[t,j,])
as <- matrix(c(0),ne,1)
for (i in 1:nm)
{
as <- as + au[t,,i,j] * (utpr[t,j,i] / pr[(t + 1),j])
}
auc[(t + 1),,j] <- as
Ps <- matrix(c(0),ne,ne)
for (i in 1:nm)
{
Ps <-
Ps + (Pu[t,,,i,j] + (auc[(t + 1),,j] - au[t,,i,j]) %*% (t(auc[(t + 1),,j] -
au[t,,i,j]))) * (utpr[t,j,i] / pr[(t + 1),j])
}
Puc[(t + 1),,,j] <- Ps
}
}
###
### If there are missing values we have a separate time loop... (for no reason)
###
if (any(is.na(y[t,])))
{
### ----------------- debug code ------------------- ###
if(debug)
cat(sprintf("missing value(s) at t = %d\n", t));
### ----------------- end debug ------------------- ###
for (j in 1:nm)
{
for (i in 1:nm)
{
# KALMAN FILTER STEP
# PREDICTION EQUATIONS
ap[t,,i,j] <- c[1,,1,j] + H[1,,,j] %*% (auc[t,,i])
Pp[t,,,i,j] <-
H[1,,,j] %*% Puc[t,,,i] %*% t(H[1,,,j]) + G[1,,,j] %*% K[1,,,j] %*% t(G[1,,,j])
# ONE-STEP-AHEAD PREDICTION ERROR CAN NOT BE CALCULATED
# BUT THE ASSOCIATED COVARIANCE MATRIX CAN
F[t,,,i,j] <- W[1,,,j] %*% Pp[t,,,i,j] %*% t(W[1,,,j]) + R[1,,,j]
# THUS, THE UPDATED STATE IS IDENTICAL TO THE PREDICTED STATE
# AND THE ASSOCIATED COVARIANCE MATRIX IS CALCULATED AS USUAL
au[t,,i,j] <- ap[t,,i,j]
Pu[t,,,i,j] <-
Pp[t,,,i,j] + Pp[t,,,i,j] %*% t(W[1,,,j]) %*% solve(F[t,,,i,j]) %*% W[1,,,j] %*% Pp[t,,,i,j]
# END OF KALMAN FILTER STEP
# BEGIN HAMILTON FILTER STEP
# CONDITIONAL TRANSITION PROBABILITY FOR S(t-1)=i to S(t)=j
tpr[t,j,i] <- p[j,i] * pr[t,i]
# CONDITIONAL DENSITY CAN NOT BE CALCULATED
}
}
# CURRENT OCCASION DOES NOT CONTRIBUTE TO THE LIKELIHOOD FUNCTION
utpr[t,,] <- tpr[t,,]
# END OF THE HAMILTON FILTER STEP
# BEGIN COLLAPSE STEP to reduce the number of posteriors
for (j in 1:nm)
{
pr[(t + 1),j] <- sum(utpr[t,j,])
as <- matrix(c(0),ne,1)
for (i in 1:nm)
{
as <- as + au[t,,i,j] * (utpr[t,j,i] / pr[(t + 1),j])
}
auc[(t + 1),,j] <- as
Ps <- matrix(c(0),ne,ne)
for (i in 1:nm)
{
Ps <-
Ps + (Pu[t,,,i,j] + (auc[(t + 1),,j] - au[t,,i,j]) %*% (t(auc[(t + 1),,j] -
au[t,,i,j]))) * (utpr[t,j,i] / pr[(t + 1),j])
}
Puc[(t + 1),,,j] <- Ps
}
}
} # END OF TIME LOOP
Lp<-if(is.nan(L)) 1e9 else min(-L,1e9)
attr(Lp,"jd") <- jd
attr(Lp,"all") <- list(jd=jd, Puc=Puc, Pu=Pu, Pp=Pp, v=v, F=F, pr=pr, tpr=tpr, utpr=utpr, au=au, ap=ap, auc=auc, md=md)
Lp
}
Try the mskf package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.