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),'\\')
}