R help - Aug 2008 - Maximum likelihood estimation

Hello, 

I am struggling for some time now to estimate AR(1) process for commodity price
time series. I did it in STATA but cannot get a result in R.

The equation I want to estimate is: p(t)=a+b*p(t-1)+error 
Using STATA I get 0.92 for a, and 0.73 for b. 

Code that I use in R is: 
p<-matrix(data$p) # price at time t
lp<-cbind(1,data$lp) # price at time t-1

                                             
mle <- function(theta) {
  sigma2<-theta[1]
  b<- theta[-1]
  n<-length(p)
  e<-p-lp%*%b
  logl<- -(n/2)*log(sigma2)-((t(e)%*%e)/(2*sigma2))
  return(-logl)
  }


out <- optim(c(0,0,0),mle, method = "L-BFGS-B",
    lower = c(0, -Inf, -Inf),
    upper = c(Inf, Inf, Inf))

The "result" I get is: " Error in optim(c(0, 0, 0), mle, method =
"L-BFGS-B", lower = c(0, -Inf,:L-BFGS-B needs finite values of
'fn'"

Can somebody spot the mistake?

Many thanks, 

Jurica Brajkovic

Jurica Brajkovi? wrote:
> Hello, 
>
> I am struggling for some time now to estimate AR(1) process for commodity
price time series. I did it in STATA but cannot get a result in R.
>
> The equation I want to estimate is: p(t)=a+b*p(t-1)+error 
> Using STATA I get 0.92 for a, and 0.73 for b. 
>
> Code that I use in R is: 
> p<-matrix(data$p) # price at time t
> lp<-cbind(1,data$lp) # price at time t-1
>
>                                              
> mle <- function(theta) {
>   sigma2<-theta[1]
>   b<- theta[-1]
>   n<-length(p)
>   e<-p-lp%*%b
>   logl<- -(n/2)*log(sigma2)-((t(e)%*%e)/(2*sigma2))
>   return(-logl)
>   }
>
>
> out <- optim(c(0,0,0),mle, method = "L-BFGS-B",
>     lower = c(0, -Inf, -Inf),
>     upper = c(Inf, Inf, Inf))
>
> The "result" I get is: " Error in optim(c(0, 0, 0), mle,
method = "L-BFGS-B", lower = c(0, -Inf,:L-BFGS-B needs finite values
of 'fn'"
>
> Can somebody spot the mistake?
>
> Many thanks, 
>
> Jurica Brajkovic 
>
>   
As far as I can see, the first element of the vector of starting values
(sigma2) is 0 and you are taking log of it....
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
>   

-- 
   O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark      Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)              FAX: (+45) 35327907

Hi R-experts,
I'm new to R in mle. I tried to do the following but just couldn't get
it
right. Hope someone can point out the mistakes. thanks a lot.

t <- c(1:90)
y <-
c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,

234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,

401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,

464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
	473,473,473,475,475,475,475)
fr <- function(a, b, p, lambda){
l <- 0.5*(lambda + b*p + (1-p)*(lambda-b))
l^2 > lambda*b*p
delta <- sqrt(abs(l^2 - b*p*lambda))
mt <- a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
logl <- sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
return(-logl)
}
optim(c(15,0.01,0.5,0.01),fr, method="L-BFGS-B", 
lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
Inf),control=list(fnscale=-1))

-- 
View this message in context:
http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19304249.html
Sent from the R help mailing list archive at Nabble.com.