[R] nlminb error
Bert Gunter
bgunter.4567 at gmail.com
Mon Apr 10 06:43:39 CEST 2017
Whether you have converged to a global optimum or not depends on the
nature of the likelihood surface. Have you tried different starting
values? As for the warning, I leave that to Prof. Nash, as he *is*
(way) more knowledgeable. However, I suspect the answer is yes, it is
concerning. Are you at a boundary?
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 Sun, Apr 9, 2017 at 5:55 PM, Rodrigo Andres Cristancho Castellanos
<cristanchocc at gmail.com> wrote:
> Hi Bert and JC, thanks for the quick response.
>
> Thinking about what Bert mentions as posible "issues" I think the three of
> them are feasible explanations of the error.
>
> Anyway, I believe I have already found a work around, I just replace the
> function "nlminb" with the function "optim". Even though, the warning
> persits the result seems to converge without other issues, given that it
> indicates that the convergence is succesfully completed.
>
> The only question that comes to me with this lucky strike, is, do I need
> keep worring about the warning ?
>
> Thank you both for your help.
>
> Regards,
>
> Rodrigo Cristancho
>
>
>
> 2017-04-09 12:18 GMT-05:00 Bert Gunter <bgunter.4567 at gmail.com>:
>>
>> 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
>> <cristanchocc at 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 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