R help - Nov 2014 - optimx

Hello All,

I need some help with optimx, mostly due to my apparent lack of imagination than
optimx itself.  Below is the pertinent code for which I have a question.  I am
fitting to the term structure of swap rates per Cox, Ingersoll, and Ross.  As
you can see the objective function is CIRTune.  I have set up the lower and
upper bounds of the model for each kappa, lambda, and theta.  Everything works
as expected.

However, there is an additional constraint.
2 * kappa * lamba >= sigma^2

My question is as follows:  should the constraint be worked into the function or
should I add an additional parameter in the function that is controlled by lower
and upper?

 #Objective function
  CIRTune <- function(param = numeric(), shortrate = numeric(), sigma = .015,
cfmatrix = matrix(), matmatrix = matrix()){
    kappa =   param[1]
    lambda = param[2]
    theta =     param[3]

    Disc <- CIRBondPrice(kappa = kappa, lambda = lambda, theta = theta,
shortrate = shortrate, T= matmatrix,  step = 0, sigma = sigma)
    
   CIRTune <- sqrt((sum(colSums((cfmatrix * Disc))^2))/ncol(matmatrix))
    return(CIRTune)
  }
  
  # Fit the model to the market   
  fit <- optimx(par = c(.1, .003, .03), 
                fn = CIRTune, 
                method = "L-BFGS-B",
                lower = c(rep(.001,3),
                upper = rep(1, 3), 
                shortrate = shortrate,
                sigma = .015,
                cfmatrix = CIR.CF.Matrix, 
                matmatrix = CIR.Mat.Matrix)  
  close(CalCIR1)
  return(fit)


Thank-you in advance,
Glenn
	[[alternative HTML version deleted]]

Glenn Schultz <glennmschultz <at> me.com> writes:
> 
> Hello All,
 
> I need some help with optimx, mostly due to my apparent lack of
> imagination than optimx itself.  Below is the pertinent code for
> which I have a question.  I am fitting to the term structure of swap
> rates per Cox, Ingersoll, and Ross.  As you can see the objective
> function is CIRTune.  I have set up the lower and upper bounds of
> the model for each kappa, lambda, and theta.  Everything works as
> expected.
 
> However, there is an additional constraint.
> 2 * kappa * lamba >= sigma^2
 
> My question is as follows: should the constraint be worked into the
> function or should I add an additional parameter in the function
> that is controlled by lower and upper?
  This is a general question about constrained optimization, not
just about optimx (e.g. it would apply to the L-BFGS-B method for
base::optim() as well).  In general the best way to handle this
(I think) would be reparameterize your model slightly, substituting
gamma = 2*kappa*lambda-sigma^2 for either lambda or kappa.  Then
you know that gamma must be >=0, and you can compute the appropriate
value of (e.g.) lambda = (gamma+sigma^2)/(2*kappa) on the fly from
the specified values of kappa and gamma.  This will automatically
satisfy the lambda>0 constraint (if gamma, sigma, kappa are all >0),
but you need to check what's necessary to satisfy the other
(lambda<1) constraint; I haven't worked all the way through that
("left as an exercise"), but if it can't definitely be satisfied
by an independent constraint on gamma then you're in trouble
(you have to work out what values of kappa, sigma, and gamma
could lead to lambda>1, then see if it's possible to prevent them
by assigning independent constraints).

If you *are* in trouble, then you need to find an optimizer
that can handle nonlinear inequality constraints.  I think the nloptr
package might have one, or you can do this by hand (1) exactly by
Lagrange multipliers or (2) approximately by imposing a nonlinear penalty
when the constraint is violated.
>  #Objective function
>   CIRTune <- function(param = numeric(), 
> shortrate = numeric(), sigma = .015, cfmatrix = matrix(),
> matmatrix = matrix()){
>     kappa =   param[1]
>     lambda = param[2]
>     theta =     param[3]
> 
>     Disc <- CIRBondPrice(kappa = kappa, lambda = lambda, 
> theta = theta, shortrate = shortrate, T= matmatrix, 
> step = 0, sigma = sigma)    
> 
>    CIRTune <- sqrt((sum(colSums((cfmatrix * Disc))^2))/ncol(matmatrix))
>     return(CIRTune)
>   }
> 
>   # Fit the model to the market   
>   fit <- optimx(par = c(.1, .003, .03), 
>                 fn = CIRTune, 
>                 method = "L-BFGS-B",
>                 lower = c(rep(.001,3),
>                 upper = rep(1, 3), 
>                 shortrate = shortrate,
>                 sigma = .015,
>                 cfmatrix = CIR.CF.Matrix, 
>                 matmatrix = CIR.Mat.Matrix)  
>   close(CalCIR1)
>   return(fit)