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.l...@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]]
______________________________________________
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.
