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