Using R to find MLE's

#1
I have a model like this (Nelson and Siegel 1987)



where tau and y are the maturities and yields, respectively, given to me in my data file..

y=
4.863,5.662,6.41,6.864,7.153,7.352,7.409,7.474,7.503,7.644,7.676,7.701,7.674,7.668,7.665,7.741,7.743,7.742
tau=
1,3,6,9,12,15,18,21,24,30,36,48,60,72,84,96,108,120

I firstly need to find the MLE of m which maximises the likelihood function and then I can easily find the b1, b2 and b3 constants using this m value via least squares estimation..

But does anybody know how I can go abouts finding the MLE of m and if you could help with providing r code for it, I would appreciate that a lot. I have been pulling my hair out for the past week now trying to do it :)

P.S. I know that m must be in between 0 and 0.4
 
Last edited:

Dason

Ambassador to the humans
#2
You need a likelihood to find an MLE. I don't see any distributional assumptions here. I also don't see what y(m) is supposed to be - are you claiming that the data is exactly determined by the other parameters? In short - I'm kind of confused as to what you're doing.
 
#3
Hi, thanks for the quick reply.

I need to find the m (lambda) value that maximizes the likelihood, where the likelihood is:

[ (n/2)*log(sigma) ] + [ sum(z^2)/2*(sigma^2) ]

where z = y - beta1 - beta2x1 - beta3x2

note that x1 and x2 depend on m (as seen in the picture above, involving the tau and exponentials)

And of course I need to do this in R, but I am struggling to find out which m value maximizes this function, where m is in between 0 and 0.4 e.g. I could do seq(0,0.4,length=50) for example, but I still don't understand how exactly I can make R give me the MLE of m.

y is the yield and depends on m, as well as b1, b2, b3 and tau, i am given the y values and tau values and need to find the MLE of m, which will help me get b1, b2, b3, using lm(y~x1+x2)
 
#5
the L has sum(z^2) inside it, z contains x1 and x2 and x1 and x2 depend on m (first picture in first post, x1 = 1-exp, etc), y is the data yes, i have 18 y values and 18 tau values

here is a little function i did:
myfun = function(U)
{
#U[1] = beta1
#U[2] = beta2
#U[3] = beta3
#U[4] = sigma
#U[5] = lambda
n = length(y) #=18
x1 = (1-(exp(-U[5]*tau)))/(U[5]*tau)
x2 = ((1-(exp(-U[5]*tau)))/(U[5]*tau)) - (exp(-U[5]*tau))
z = y - U[1] - U[2]*x1 - U[3]*x2
L = (n/2)*log(U[4]) + (sum(z^2)/(2*(U[4]^2)))
L
}

m = lambda by the way
 
Last edited:

BGM

TS Contributor
#6
I see you have the normal assumption on the yield curve; in this case you are assuming the Nelson Siegel model the mean of the yield curve, and so you can just simultaneously find out the MLE of all parameters (Note that the maximum likelihood estimator of the coefficients of linear regressors in the mean of normal likelihood is equivalent to the least-square estimate as you have the square inside the exponent)

Or equivalently the least-square estimate can be expressed as a function of \( m \) in closed form, in which you can plug in and maximize in \( m \).
 
#8
Im sorry guys but I have no idea how to do this with R (finding mle of m) so I ask if you can spend a little bit of your time to try it with R.

Regards
 

bryangoodrich

Probably A Mammal
#9
First, when writing code, enclose it in [noparse]
Code:
 ... your code here ...
[/noparse] BB tags so we can see it appropriately as

Code:
 ... your code here ...
Just looking over your function, you haven't defined tau anywhere that I can see. Either you're referencing it outside of the function (bad programming) or you forgot to specify/pass it into the function.
 
#10
tau = 1,3,6,9,12,15,18,21,24,30,36,48,60,72,84,96,108,120, already done this.. but i just don't know how to find the mle of m, where on earth do i find the m? after doing that function i did above, i have no clue what to do next
 
#11
For example, say if I did:

Code:
> m
 [1] 0.01000000 0.01795918 0.02591837 0.03387755 0.04183673 0.04979592 0.05775510 0.06571429 0.07367347 0.08163265 0.08959184
[12] 0.09755102 0.10551020 0.11346939 0.12142857 0.12938776 0.13734694 0.14530612 0.15326531 0.16122449 0.16918367 0.17714286
[23] 0.18510204 0.19306122 0.20102041 0.20897959 0.21693878 0.22489796 0.23285714 0.24081633 0.24877551 0.25673469 0.26469388
[34] 0.27265306 0.28061224 0.28857143 0.29653061 0.30448980 0.31244898 0.32040816 0.32836735 0.33632653 0.34428571 0.35224490
[45] 0.36020408 0.36816327 0.37612245 0.38408163 0.39204082 0.40000000
And these are my possible values of m,

How would I link this to the function above? Wouldn't I need to create a loop?
 

bryangoodrich

Probably A Mammal
#12
If your function is step up correctly, no. You should have no reason for a loop. Consider minimizing a quadratic.

Code:
x <- seq(-10, 10)
f <- function(x) x^2 + 2*x - 1
min(f(x))  # = -2
If all you're doing is looking for the maximum value returned by your function for each m, then your function ought to return a similar vector of values for each m. This is called vectorization. It is implicit in the example above for how the min function in R works. The same applies to sum and other arithmetic. It evaluates a vector. Your function should be designed in similar thought: to evaluate a vector. Frankly, I see 3 different functions within your function. You should consider breaking it down:

Code:
myfun <- function(U) {
  n  <- length(y) #=18
  x1 <- function(lam, tau) (1 - exp(-lam*tau)) / (lam*tau)
  x2 <- function(lam, tau) ((1 - exp(-lam*tau)) / (lam*tau)) - exp(-lam*tau)
  z  <- y - U[1] - U[2]*x1(U[5], tau) - U[3]*x2(U[5], tau)
  L  <- (n/2) * log(U[4]) + (sum(z^2) / (2*U[4]^2))
  return(L)
}
But as I already said, you haven't passed tau or y into this function. I can only assume they are defined outside of the function. While R will make the appropriate reference to those vectors, your function is ill-written as it stands. You also probably shouldn't pass it one object that contains several. Pass each thing individually, or at least group them logically, like

Code:
myfun <- function(m, y, tau, betas, sigma) {
  n  <- length(y) #=18
  x1 <- function(lam, tau) (1 - exp(-lam*tau)) / (lam*tau)
  x2 <- function(lam, tau) ((1 - exp(-lam*tau)) / (lam*tau)) - exp(-lam*tau)
  z  <- y - betas[1] - betas[2]*x1(m, tau) - betas[3]*x2(m, tau)
  L  <- (n/2) * log(sigma) + (sum(z^2) / (2*sigma^2))
  return(L)
}
This, to me, appears much more readable, and you only need to pass it a grouped element of the coefficients (betas).
 

bryangoodrich

Probably A Mammal
#14
If you want to use that function to linearly optimize it, you either have to know how to code your own linear optimization function, do something analytic, or hope a pre-defined package out there will work with your function. I know a bit about numerical analysis, but I've avoided linear programming classes so far, so I don't know the algorithms they use to do it. My guess is that your problem has a more analytical approach, since you are just fitting a linear model: \(Y = \beta_0 + \beta_1 X^* + \beta_2 X^{**}\), where each of those X's are themselves functions of m (tau is given, so it appears exogenous to me). I'm too rusty on this stuff to be much help at the moment.

But based on what I said before, if everything but m is determined as inputs, then you need to write your function so you can say "f(m)" where m is a vector of values (e.g., m = seq(0, 0.4, by = 0.001)), and get back a vector of values. Then you simply have to pick out the maximum, which is exactly what I showed you in my prior example with optimizing the quadratic. If you don't have everything as inputs, the problem is obviously more difficult. Either linearly optimize the function as is, or figure out how to get the rest of the inputs based on what you know.

Though, if you want something R has for dealing with MLE's, there's the mle function in the stat4 package, and the nlm function in the nlm package. It's incumbent upon you to figure out how to put your problem into a way those functions can operate on it, but my guess is that there is a much easier solution given the information you have.
 
#15
All I want to do is find the m that maximizes my L :(

I am pulling my hair out seriously, I just don't get this at all!

I have all my possible values for m, 50 of them, I just want to know which one of them maximizes the likelihood function and I have absolutely no clue how to do it in R, nothing it making sense to me.

Sure, I can do:

myfun(1,2,3,4,5) where myfun(beta1,beta2,beta3,sigma,m) and it gives me an answer for the likelihood L, but I can't do that, I need to find the exact value of m from these values:

Code:
[1] 0.01000000 0.01795918 0.02591837 0.03387755 0.04183673 0.04979592 0.05775510 0.06571429 0.07367347 0.08163265 0.08959184 0.09755102
[13] 0.10551020 0.11346939 0.12142857 0.12938776 0.13734694 0.14530612 0.15326531 0.16122449 0.16918367 0.17714286 0.18510204 0.19306122
[25] 0.20102041 0.20897959 0.21693878 0.22489796 0.23285714 0.24081633 0.24877551 0.25673469 0.26469388 0.27265306 0.28061224 0.28857143
[37] 0.29653061 0.30448980 0.31244898 0.32040816 0.32836735 0.33632653 0.34428571 0.35224490 0.36020408 0.36816327 0.37612245 0.38408163
[49] 0.39204082 0.40000000
Which will maximize the likelihood function.

I am asking please if somebody has the time and patience, try to help me with do this question. I will even pay you if you want!

Nobody in my class knows what to do, they are all struggling with finding the m value that maximizes L
 
Last edited:

TheEcologist

Global Moderator
#16
Will it help if you have a simple example?

Lets say counts of malaria mosquitoes close to a stagnant pond follows this function (regardless of direction); [TEX]\tau \lambda e^{-\lambda r}[/TEX]. Where r is distance in meters, [TEX]\lambda[/TEX] is the decay of biting mosquitoes with distance and [TEX]\tau[/TEX] is the total amount of dispersed female mosquitoes. At each point the number of biting mosquitoes is Poisson distributed around an expected value given by the function above.

I want to find the value of [TEX]\lambda[/TEX] and [TEX]\tau[/TEX] for a given pond from which I have taken samples at different distances. The r code would look like this;


Code:
# first generate some simulated data with 
lambda = 1/50; tau =45000
#sample distances
r=1:100
# measurements
expectedbites=tau*lambda*exp(-lambda*r)
#simulated actual bites
bites=rpois(100,expectedbites)

# now lets back fit the values
# our log-likelihood function for parameters p is
LL=function(p) {
predictions=p[1]*p[2]*exp(-p[2]*r)
sum(-dpois(bites,predictions,log=T))
}

#select some start values
startvals=c(10000,1/20) 
# now optimize our LL function
optim(startvals, LL)
And you will see that we have found pretty good estimates for this problem!

Hope this helps you to adapt your code.

TE
 

BGM

TS Contributor
#17
Yes the function optim is frequently use in R for numerical optimization.

I think the idea from OP should refer to the bisection method, and actually it is not that hard to code that numerical method by yourself also. Of course you have a nice built-in function optim to do it for you :)
 

TheEcologist

Global Moderator
#20
How do I know what start values to use for b1, b2, b3, sigma and lambda? Because if I choose different ones, I will get different answers
That sounds like an over-parameterization problem, with many local (likelihood) maxima [or local minima for the log-likelihood]. Best thing to do is to try a large range of different start values until you are sure you can find the solution that minimizes the log-likelihood.

You can also try if a simpler model doesn't do a better job.