Hi,
I'm trying to estimate the parameters of a state space model of the
following form
measurement eq:
z_t = a + b*y_t + eps_t
transition eq
y_t+h = (I -exp(-hL))theta + exp(-hL)y_t+ eta_{t+h}.
The problem is that the distribution of the innovations of the transition
equation depend on the previous value of the state variable.
To be exact: y_t|y_{t-1} ~N(mu, Q_t) where Q is a diagonal matrix with
elements equal to
Q_{i,t} =
sigma_i*(1-exp(-kappa_i*h)/kappa_i*(theta_i/2*(1-exp(kappa_i*h)+exp(-kappa_i*h)y_{t-1,i}
The fkf returns the filtered states variables so y_{t-1,i} is available. I
just can't figure out how to write my program in such a way that this
information is included and updated in the state space model for each
iteration in optim. Any suggestions on how to solve this are much
appreciated.
Thank you,
Kristian.
Below my program where the variance matrices are just identity matrices and
the data is just random numbers. I used the example from the FKF package as
framework.
library(FKF) #loading Fast Kalman Filter package
library(Matrix) # matrix exponential package
library(BB)
library(alabama)
x <- matrix(abs(rnorm(400)), nrow=10, ncol=40)
m <- 2 # m is the number of state variables
n <- ncol(x) # is the length of the observed sample
d <- nrow(x) # is the number of observed variables.
h <- 1/52
## creating state space representation of 2-factor CIR model
CIR2ss <- function(K_1, K_2, sigma_1, sigma_2, lambda_1, lambda_2, theta_1,
theta_2, delta_0, delta_1, delta_2) {
## defining auxilary parameters,
phi_11 <- sqrt((K_1+lambda_1)^2+2*sigma_1^2*delta_1)
phi_21 <- sqrt((K_2+lambda_2)^2+2*sigma_2^2*delta_2)
phi_12 <- K_1+lambda_1+phi_11
phi_22 <- K_2+lambda_2+phi_21
phi_13 <- 2*K_1*theta_1/sigma_1^2
phi_23 <- 2*K_2*theta_2/sigma_2^2
a <- array(1, c(d,n))
phi_14 <- numeric(d)
phi_24 <- numeric(d)
b1 <- numeric(d)
b2 <- numeric(d)
for(t in 1:d){
phi_14[t] <- 2*phi_11+phi_12*(exp(phi_11*t)-1)
phi_24[t] <- 2*phi_21+phi_22*(exp(phi_21*t)-1)
a[t] <- delta_0 - phi_13/(t)*log(2*phi_11*exp(phi_12*(t)/2)/phi_14[t])-
phi_23/(t)*log(2*phi_21*exp(phi_22*(t)/2)/phi_24[t])
b1[t] <- (delta_1/(t))*2*(exp(phi_11*(t))-1)/phi_14[t]
b2[t] <- (delta_2/(t))*2*(exp(phi_21*(t))-1)/phi_24[t]
}
b <- array(c(b1,b2), c(d,m,n))
j <- -matrix(c(K_1, 0, 0, K_2), c(2,2))*h
explh <- as.matrix(expm(j))
Tt <- array(explh, c(m,m,n)) #array giving the factor of the transition
equation
Zt <- b #array giving the factor of the measurement equation
ct <- a #matrix giving the intercept of the measurement equation
dt <- as.matrix((diag(m)-explh)%*%c(theta_1, theta_2)) #matrix giving
the intercept of the transition equation
GGt <- array(diag(d), c(d,d,1)) #array giving the variance of the
disturbances of the measurement equation
HHt <- diag(m) #array giving the variance of the innovations of the
transition equation
a0 <- c(0.5, 0.5) #vector giving the initial value/estimation of the
state variable
P0 <- matrix(1e6, nrow = 2, ncol = 2) # matrix giving the variance of a0
return(list(a0 = a0, P0 = P0, ct = ct, dt = dt, Zt = Zt, Tt = Tt, GGt =
GGt,
HHt = HHt))
}
## Objective function passed to optim
objective <- function(parm, yt) {
sp <- CIR2ss(parm["K_1"], parm["K_2"], parm["sigma_1"], parm["sigma_2"],
parm["lambda_1"], parm["lambda_2"],
parm["theta_1"], parm["theta_2"], parm["delta_0"],
parm["delta_1"], parm["delta_2"])
ans <- fkf(a0 = sp$a0, P0 = sp$P0, dt = sp$dt, ct = sp$ct, Tt = sp$Tt,
Zt = sp$Zt, HHt = sp$HHt, GGt = sp$GGt, yt = yt)
return(-ans$logLik)
}
parm <- c(K_1 = 0.6048, K_2 = 0.656, sigma_1 =0.1855, sigma_2 =0.5524,
lambda_1 =0.55, lambda_2 =0.009, theta_1 = 0.173, theta_2 = 0.12,
delta_0 =0.686, delta_1 =-0.003, delta_2=0.0025)
hin <- function(parm, yt){
h<- numeric(2)
h[1] <- 2*parm[1]*parm[7]-parm[3]^2
h[2] <- 2*parm[2]*parm[8]-parm[4]^2
h
}
##optimizing objective function
fit <- auglag(par = parm, fn = objective, yt = x, control.outer =
list(method = "CG"), hin = hin)
print(fit)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
