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.