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]]
I found an old thread on R-Sig-Finance with the same problem and a possible solution https://stat.ethz.ch/pipermail/r-sig-finance/2007q2/001362.html I've used with success a few times but it seems a bit slow. If anyone has a better way of modelling state-dependent volatility using one of the available packages please let me know. Kristian 2011/11/12 Kristian Lind <kristian.langgaard.lind@gmail.com>> 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]]