SCM

Forum: help

Monitor Forum | Start New Thread Start New Thread
RE: 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

Thanks to:
Vienna University of Economics and Business Powered By FusionForge