library(XNomial) #Estimating base-matrix -- STARTING HERE P_hat=matrix(0,nrow=N,ncol=N) #sum up all contingency matrices, because of time-homogenuity sumcontin = Reduce('+',contin) for(i in 1:N){ for(j in 1:N){ P_hat[i,j] = sumcontin[i,j]/sum(sumcontin[i,]) } } basematrix-P_hat #Estimating base-matrix -- ENDING HERE #Goodness of fit test, Pearson Chi Square for the base matrix X2 = 0 for (k in 1:(T-1)){ for (i in 1:N){ Oi = sum(contin[[k]][,i]) u=vector() for (j in 1:N){ u = c(u,sum(contin[[k]][j,])) } Ei = (u%*%P_hat)[i] X2 = X2 + (Oi-Ei)^2/Ei } } X2 1-pchisq(X2, N*(T-N)) #Estimate transition probabilities for single firm -- START HERE #Finding contingency table of single firm f f=12 continsf = matrix(0,nrow=N,ncol=N) for(k in 2:T){ continsf[states[f,k-1],states[f,k]] = continsf[states[f,k-1],states[f,k]]+1 } continsf #Use maximum likelihood estimation P_hatsf = matrix(0,nrow=N,ncol=N) for(i in 1:N){ for(j in 1:N){ P_hatsf[i,j] = continsf[i,j]/sum(continsf[i,]) } } #Estimate transition probabilities for single firm -- ENDING HERE #Goodness of fit test (Multinomial) for single firm with base matrix estimation M=vector() for(i in 1:N){ observed = continsf[i,] if(sum(observed)>0){ expected = P_hat[i,] V=N j=1 while(j<=V){ if(expected[j]==0 && observed[j]==0){ observed=observed[-j] expected=expected[-j] j=j-1 V=V-1 } j=j+1 } mul = xmulti(observed,expected,detail=2); M = c(M,mul\$pProb) } } M sum(M)/length(M) #Improve goodness of fit by linear combination approach beta = seq(0.01,0.2,by=0.01) P_hat_refined = list() for(k in 1:(length(beta))){ P_hat_refined[[k]] = P_hat+beta[k]*P_hatsf normal = vector() for(i in 1:N){ normal[i]=sum(P_hat_refined[[k]][i,]) } for(i in 1:N){ P_hat_refined[[k]][i,]=P_hat_refined[[k]][i,]*1/normal[i] } } #Test goodness of fit for the new refined P_hat estimation k_test = 2 M=matrix(0,nrow=length(beta),ncol=N) sink("output.txt") for(k in 1:length(beta)){ for(i in 1:N){ observed = continsf[i,] if(sum(observed)>0){ expected = P_hat_refined[[k]][i,] V=N j=1 while(j<=V){ if(expected[j]==0 && observed[j]==0){ observed=observed[-j] expected=expected[-j] j=j-1 V=V-1 } j=j+1 } mul = xmulti(observed,expected,detail=2); M[k,i]=mul\$pProb } } } sink() M M_avg = vector() for(k in 1:length(beta)){ M_avg[k]=sum(M[k,])/length(M[k,]) } latex = matrix(0,nrow=length(beta),ncol=N+2) for(k in 1:length(beta)){ latex[k,1]=paste(beta[k],' &') for (i in 2:(N+1)){ latex[k,i]=paste(toString(round(M[k,i-1],4)),' &') } latex[k,N+2]=paste(toString(round(M_avg[k],4)),' \\') } latex_P = matrix(0,nrow=N,ncol=N) for(i in 1:N){ latex_P[i,1]=paste(i,'&',round(P_hat_refined[[10]][i,1],4),'&') for (j in 2:(N-1)){ latex_P[i,j]=paste(round(P_hat_refined[[10]][i,j],4),'&') } latex_P[i,N]=paste(round(P_hat_refined[[10]][i,N],4),'\\') }