I suspect, as you hinted, there's little to no hope that anyone will
be willing or able to navigate your code. **Usually** (whatever that
means!) these sorts of problems can be traced back to
overparameterization -- too few data, which could also mean a lot of
"correlated" data, chasing too many parameters (leading to ridges in
the surface or problems near the boundaries). This is the bane of
numerical optimizers.
Try fitting simpler models with fewer parameters if possible, at least
to see if convergence can be achieved. Better starting values, other
optimizers, and/or use of analytical gradients and hessians can also
sometimes help.
Sorry I can't be more specific. Perhaps someone more knowledgeable
than I can be.
Cheers,
Bert
Bert Gunter
"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
On Sat, Apr 8, 2017 at 5:21 PM, Rodrigo Andres Cristancho Castellanos
<cristanch...@gmail.com> 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@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.
______________________________________________
R-help@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.
