Forum: help
Monitor Forum | Start New ThreadRE: need some help about how to use maxLik(2) [ Reply ] By: Arne Henningsen on 2011-03-30 17:21 | [forum:4209] |
What do you mean with "attach" a file? Please ask questions that are not specific to maxLik (e.g. reading data files) on the R help mailing list. |
need some help about how to use maxLik(2) [ Reply ] By: Jessie Xu on 2011-03-30 15:35 | [forum:4207] |
I don't know why i can not attach files. So my code is here: # Test Case1: # Y(i)=b(1)*Y(i-1) # There is only 3 states of price change ticks, -1,0,1 # sigma=1 # This code is going to estimate para={beta1,alph1,alph2} library(maxLik) sample<-as.matrix(read.table("sample1.txt",sep=","),,byrow=TRUE) n.sample<-dim(sample)[1] # number of sample n.indX<-dim(sample)[2] # number of independent variables state<-c(-1,0,1) # state for number ticks we are counting for distribution n.states<-length(state) indexMatrix<-matrix(0,n.sample,2) YMatrix<-matrix(0,n.sample,length(state)) for(i in 1:n.sample) #need to simplify by matrix function { if(sample[i,1]==-1) indexMatrix[i,]<-c(i,1) else if(sample[i,1]==0) indexMatrix[i,]<-c(i,2) else indexMatrix[i,]<-c(i,3) } YMatrix[indexMatrix]<-1 sigma=1 loglik<-function(para){ temp1=(para[2]-para[1:1]%*%t(sample))/sigma temp2=(para[3]-para[1:1]%*%t(sample))/sigma sum(YMatrix[,1]*log(pnorm(temp1)) +YMatrix[,2]*log(pnorm(temp2)-pnorm(temp1)) +YMatrix[,3]*log(1-pnorm(temp2))) value=sum(YMatrix[,1]*log(pnorm(temp1)) +YMatrix[,2]*log(pnorm(temp2)-pnorm(temp1)) +YMatrix[,3]*log(1-pnorm(temp2))) print(para) print(value) } grad<-function(para) { g <- matrix(0,n.sample,length(para)) temp1=(para[2]-para[1:1]%*%t(sample))/sigma temp2=(para[3]-para[1:1]%*%t(sample))/sigma g[,1]=t(-sample/sigma)*(-1*YMatrix[,1]*dnorm(temp1)/pnorm(temp1) -1*YMatrix[,2]*(dnorm(temp2)-dnorm(temp1))/(pnorm(temp2)-pnorm(temp1)) +YMatrix[,3]*dnorm(temp2)/(1-pnorm(temp2))) g[,2]=t((1/sigma)*(YMatrix[,1]*dnorm(temp1)-1*YMatrix[,2]*dnorm(temp1)/(pnorm(temp2)-pnorm(temp1)))) g[,3]=t((1/sigma)*(YMatrix[,2]*dnorm(temp2)/(pnorm(temp2)-pnorm(temp1))-1*YMatrix[,3]*dnorm(temp2)/(1-pnorm(temp2)))) return(numericGradient(loglik,c(0.1,-2,1))) } #output1<-nlm(loglik,para<-c(0.1,-1,1),hessian=TRUE) output <- maxBHHH(loglik,grad,, para<-c(1,-2,1)) #output3<-maxLik(loglik,,,para<-c(0.8,-1,1),"BFGS") # works Sample1.txt -1 1 0 0 1 -1 0 1 -1 0 0 1 0 1 -1 0 0 -1 0 0 |