Problem solved.

It had something to do with calling expm(-array(c(K_1, 0, 0, K_2),
c(2,2))*h).

2011/9/22 Kristian Lind <kristian.langgaard.lind@gmail.com>

> Dear R users,
>
> When running the program below I receive the following error message:
> fit <- optim(parm, objective, yt = tyield, hessian = TRUE)
> Error in as.vector(data) :
>   no method for coercing this S4 class to a vector
>
> I can't figure out what the problem is exactly. I imagine that it has
> something to do with "tyield" being a matrix.  Any help on explaining what's
> going on and how to solve this is much appreciated.
>
> Thank you,
>
> Kristian
>
> library(FKF) #loading Fast Kalman Filter package
> library(Matrix) # matrix exponential package
>
> K_1 = 0.1156
> K_2 = 0.17
> sigma_1 = 0.1896
> sigma_2 = 0.2156
> lambda_1 = 0
> lambda_2 = -0.5316
> theta_1 = 0.1513
> theta_2 = 0.2055
>
> #test data
> tyield <- matrix(data = rnorm(200), nrow =2, ncol =100)
>
> # defining dimensions
> m <- 2 # m is the number of state variables
> n <- 100 # is the length of the observed sample
> d <- 2 # is the number of observed variables.
> theta <- c(theta_1, theta_2)
> h <- t <- 1/52 # time between observations
>
> ## creating state space representation of 2-factor CIR model follwing
> Driessen and Geyer et al.
> CIR2ss <- function(K_1, K_2, sigma_1, sigma_2, lambda_1, lambda_2, theta_1,
> theta_2)    {
>       ## defining auxilary parameters
>     phi_11 <- sqrt((K_1+lambda_1)^2+2*sigma_1^2)
>       phi_21 <- sqrt((K_1+lambda_2)^2+2*sigma_2^2)
>         phi_12 <- K_1+lambda_1+phi_11
>         phi_22 <- K_2+lambda_2+phi_12
>         phi_13 <- -2*K_1*theta_1/sigma_1^2
>         phi_23 <- -2*K_2*theta_2/sigma_2^2
>         phi_14 <- 2*phi_11+phi_21*(exp(phi_11*t)-1)
>         phi_24 <- 2*phi_12+phi_22*(exp(phi_12*t)-1)
>     phi <- array(c(phi_11, phi_21, phi_12, phi_22, phi_13, phi_23, phi_14,
> phi_24), c(4,2))
>     a <- array(0, c(d,n))
>     for(t in n:1){
>       a[,n-(t+1)] <-
> -phi_13/(n-(t+1))*log(2*phi_11*exp(phi_12*(n-(t+1))/2)/phi_14)-phi_23/(n-(t+1))*log(2*phi_21*exp(phi_22*(n-(t+1))/2)/phi_24)
>         }
>     b <- array(c(1,0,0,1,0), c(d,m,n))
>     j <- -array(c(K_1, 0, 0, K_2), c(2,2))*h
>     explh <- 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 <- (diag(m)-expm(-array(c(K_1, 0, 0, K_2), c(2,2))*h))*theta #matrix
> giving the intercept of the transition equation
>     GGt <- array(c(1,0,0,1), c(d,d,n)) #array giving the variance of the
> disturbances of the measurement equation
>     HHt <- array(c(1,0,0,1), c(m,m,n)) #array giving the variance of the
> innovations of the transition equation
>     a0 <- c(0, 0) #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"])
>   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.1156, K_2 = 0.17, sigma_1 = 0.1896, sigma_2 = 0.2156,
>           lambda_1 = 0, lambda_2 = -0.5316, theta_1 = 0.1513, theta_2 =
> 0.2055) # initial parameters
>
> ##optimizing objective function
>  fit <- optim(parm, objective, yt = tyield, hessian = TRUE)
>  print(fit)
>

	[[alternative HTML version deleted]]

