[R] nlminb error
J C Nash
profjcnash at gmail.com
Sun Apr 9 22:14:15 CEST 2017
Bert has already suggested that your code is too spaghetti to let us follow it.
I tried and got a msg that PKF isn't available for R 3.3.3.
Do you actually have a test that the function can be evaluated? That's often the main issue.
And it does help to try a few parameter sets to see if you have a reasonable function.
JN
On 2017-04-08 08:21 PM, Rodrigo Andres Cristancho Castellanos wrote:
> Hi, I´m having troubel using nlminb this is the warning that shows up.
>
> Warning: Cholesky factorization 'dpotrf' exited with status: 1
> Variance of the prediction error can not be computed.
> Warning: Cholesky factorization 'dpotrf' exited with status: 1
> Determinant of the variance of the prediction error can not be computed.
>
> I´m trying to optimize a likelihood function. The code is a little messy,
> you can find it at the end of the email.
>
> I´ll appreciate any kind of help.
>
> Regards.
>
> Rodrigo Cristancho.
>
>
>
> library(foreach)
> library(FKF)
> library(ggplot2)
> library(gridExtra)
>
> # Model setup
> #3 Factor Independent Model
>
> delta1<- c(0.006, 0.003, 0.007)
> kappa <- c(0.001, 0.002, 0.004)
> sigma <- c(0.002, 0.005, 0.003)
> rc <- c(0.00002)
> r1 <- c(0.00002)
> r2 <- c(0.5)
>
> #All other parameters
> T <- 110 # years
> delta_t <- 1 # monthly zeros observations
> dt <- 1 # monthly r simulations
> n <- T/dt # number of r simulations
> r_0 <- delta1
> measurement_error <- 0.001 # for zero-coupon rates
> m <- length(delta1) # dimension of state variables
> # maturity <- c(1/12 ,1/4, 1/2, 10) # zeros for 1 factor simulation
> # maturity <- c(1/12 ,1/4, 1/2, 1, 2, 3, 5, 7, 10) # zeros for 2 factor
> simulation
>
> maturity <- 61-xc #
> d <- length(maturity) # dimension of observations
> N <- T/delta_t # number of observations
> premubar1=list(premubar)
>
> # PARAMETER ESTIMATION
> # ----------------------------------------------------------
> # Random initialization of parameters
> init_params_if <- function()
> {
> delta1_init <<- runif(m, min=0.0, max=0.05)
> kappa_init <<- runif(m, min=0.0, max=0.05)
> sigma_init <<- runif(m, min=0.0, max=0.05)
> rc_init <<- runif(1, min=0.0, max=0.001)
> r1_init <<- runif(1, min=0.0, max=0.001)
> r2_init <<- runif(1, min=0.0, max=0.001)
> }
> # optimization parameter bounds
>
> upper_bound <- c(rep(c(1.0, 1.0, 1.0, 1.0),each=m), rep(0.1,d))
> lower_bound <- c(rep(c(0.0001, 0.0001, 0.0001,-1.0 ),each=m), rep(0.0001,d))
> actual_param <- c(kappa=kappa, delta1=delta1, sigma=sigma, rc=rep(rc,1),
> r1=rep(r1,1), r2=rep(r2,1))
>
> #
> #------------------------------------------------------------------------------------------------------------
> #
>
>
> #Kalman Filter for 3 Factor Indipendent Model
> Ind_Fact_KF <- function(delta1, kappa, sigma, rc, r1, r2, observations)
> {
> # initial state variable (a0: m x 1)
> #r_init <- as.vector(delta1) # unconditional mean of
> state variable
>
> r_init <- as.vector(c(rep(0,m)))
>
> # variance of state variable (P0: m x m)
> #P_init <- (sigma^2/(2*(delta1^3)))*diag(1,m,m) # unconditional
> variance of state variable
>
> P_init <- (sigma^2/(2*(delta1)))*diag(1,m,m)
>
> # intercept of state transition equation (dt: m x 1)
> C <- matrix(0,nrow = m,ncol = 1)
>
> # factor of transition equation (Tt: m x m x 1)
> F_ <- array(exp(-kappa*delta_t)*diag(m),dim=c(m,m,1))
>
> # factor of measurement equation (Zt: d x m x 1)
> B <-
> array(1/matrix(rep(delta1,d),d,byrow=TRUE)*(1-exp(-matrix(rep(delta1,d),d,byrow=TRUE)
> * matrix(rep(maturity,m),d))),dim=c(d,m,1))
>
> # intercept of measurement equation (ct: d x 1)
>
> A <-
> (1/2)*t(sigma^2/delta1^3)%*%t(((1/2)*(1-exp(-2*(matrix(rep(delta1,d),d,byrow=TRUE)*matrix(rep(maturity,m),d)))))
>
> -2*(1-exp(-2*(matrix(rep(delta1,d),d,byrow=TRUE)*matrix(rep(maturity,m),d))))
> -(matrix(rep(delta1,d),d,byrow=TRUE)*matrix(rep(maturity,m),d)))
>
> A <- matrix(-t(A)/maturity,nrow=length(maturity))
>
> B <- array(B[,,1]/matrix(rep(maturity,m),d),dim=c(d,m,1))
>
> # variance of innovations of transition (HHt: m x m x 1)
> Q <-
> array(sigma^2/(2*kappa)*(1-exp(-2*kappa*delta_t))*diag(m),dim=c(m,m,1))
>
> # variance of measurement error (GGt: d x d x 1)
>
> R <- array(diag(d)*(rc+r1*exp(r2)),dim=c(d,d,1))
> #R <- array(diag(d)*(rc),dim=c(d,d,1))
>
> ##Funcion de filtro de Kalman rapido con base al desarrollo del modelo
> arriba
> filtered_process <- fkf(a0=r_init, P0=P_init, dt=C, ct=A, Tt=F_, Zt=B,
> HHt=Q, GGt=R, yt=observations)
> return(filtered_process)
> }
>
> #aaaa=fkf(a0=r_i(nit, P0=P_init, dt=C, ct=A, Tt=F_, Zt=B, HHt=Q, GGt=R,
> yt=t(premubar))
> aaaa=Ind_Fact_KF(delta1, kappa, sigma, rc, r1, r2,observations=t(premubar))
> aaaa$logLik
> aaaa$Ft
>
> # Retrieve short rates using Kalman Filter
> retrieve_short_rates_if <- function(rates, optim_controls,
> lower_bound=NULL, upper_bound=NULL)
> {
> observations <- rates
> init_params_if()
> initial_param <<- c(delta1=delta1_init,kappa=kappa_init,
> sigma=sigma_init, rc=rc_init, r1=r1_init, r2=r2_init)
>
> if_KF_loglik <- function(x)
> {
> delta1 <- x[1:m]; kappa <- x[(m+1):(2*m)]; sigma <- x[(2*m+1):(3*m)];
> rc <- x[(3*m+1):(3*m+1)]; r1 <- x[(3*m+2):(3*m+2)];r2 <-
> x[(3*m+3):length(x)]
> return(-Ind_Fact_KF(delta1,kappa,sigma,rc,r1,r2,observations)$logLik)
> }
>
> # optimization of log likelihood function
> fitted_model <- nlminb(initial_param, if_KF_loglik,
> control=optim_controls, lower=lower_bound, upper=upper_bound)
> return(fitted_model)
> }
>
> rrif=retrieve_short_rates_if(rates=t(premubar),optim_controls=optim_controls,upper=upper_bound,lower=lower_bound)
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>
More information about the R-help
mailing list