dat<- read.table ("Contracts_Red.txt",header=TRUE) Anzahl<- c(4,3,3,3,3,2,2) Zielstates<- read.table("Zielstates.txt", header=FALSE) dist <- function (x) { # =1-distr. function[x1] shape<- max(0.00001,x[2]) scale<- max(0.00001,x[3]) w<- pgamma(x[1]+0.5, shape, scale, lower.tail = FALSE, log.p = FALSE) w } dist1 <- function (x) { shape<- max(0.00001,x[2]) scale<- max(0.00001,x[3]) if (x[1]<=0.5) w<- 1 else w<- pgamma(x[1]-0.5, shape, scale, lower.tail = FALSE, log.p = FALSE) w } dens <- function (x) w<- dist1(x)-dist(x) L<- function(z) { #z[1]=state, z[i] contain parameters for the transition to other states LL<- 0 for (i in 1:5199) if (dat[i,25]==z[1]){ for (m in 1:Anzahl[z[1]+1]) { x<- c(dat[i,20],z[2*m],z[2*m+1]) if (Zielstates[z[1]+1,m]==dat[i,24]) LL<- LL-log(dens(x)) else LL<- LL-log(dist(x))} if (dat[i,22]>0) for (m in 1:Anzahl[z[1]+1]){ x<- c(dat[i,22],z[2*m],z[2*m+1]) LL<- LL-log(dist(x))} } return(LL) } B <- array(0,c(7,9)) # estimated parameters and loglikelihood scores C <- array(0,c(7,8)) # standard errors of stimated parameters for (n in 1:7){ # estimation State<- n-1; Min<- 100000; Lmin<- function(y) { z<- c(State,y) r<- L(z) return(r) } dim<- 2*Anzahl[n] a<- array(0,c(dim)) anf<- array(0,c(dim)) u<- array(1,c(dim)) if (n==1) anf<- c(0.2710769, -1.028765133, 0.2339671, -8.246735903, 0.2480752, -3.103583195, 0.2623018, -8.294532) if (n==2) anf<- c(1.1270962, 0.012623693, 1.0651211, 0.013833028, 0.6512206, 0.006126045) if (n==3) anf<- c(0.9473864, 0.005466756, 1.0234273, 0.012416461, 0.6106675, 0.006318800) if (n==4) anf<- c(0.9699663, 0.006879077, 1.0550613, 0.014891264, 0.6506892, 0.005818834 ) if (n==5) anf<- c(0.8370977, 0.004940937, 0.8603592, 0.005361572, 0.5756741, 0.003256479) if (n==6) anf<- c( 0.9811662, 0.008860099, 0.9544669, 0.011138481) if (n==7) anf<- c(0.8614765, 0.005157315, 0.9324809, 0.015589898) for (iter in 1:2) {# more iterations are possible resL <- optim(par=anf,fn=Lmin,gr=NULL,method=("BFGS"),control=list(trace=1,fnscale=1,parscale=u,alpha=1,beta=0.5,gamma=2,abstol=0.00001,reltol=0.00001), hessian=FALSE) anf<- resL$par for (k in 1:(2*Anzahl[n])) if (anf[k]<0) anf[k]<- 0.01 res <- optim(par=anf,fn=Lmin,gr=NULL,method=("Nelder-Mead"),control=list(trace=1,fnscale=1,parscale=u,alpha=1,beta=0.5,gamma=2,abstol=0.00001,reltol=0.00001), hessian=FALSE) anf<- res$par for (k in 1:(2*Anzahl[n])) if (anf[k]<0) anf[k]<- 0.1 } resL <- optim(par=anf,fn=Lmin,gr=NULL,method=("BFGS"),control=list(trace=1,fnscale=1,parscale=u,alpha=1,beta=0.5,gamma=2,abstol=0.00001,reltol=0.00001), hessian=TRUE) d<-det(100*resL$hessian) if (abs(d)>10^-20) a<-diag(solve(resL$hessian)) z<- resL$par for (m in 1:(2*Anzahl[n])){ B[n,m]<- z[m] C[n,m]<- if (a[m]<0) -1 else sqrt(a[m])*10 } B[n,9]<- resL$val } print(n) print(B) print(C) # source("Gamma.R")