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 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]]