library(expm)
library(pracma)


#Simulating 200 30-state markov chain over T=30 transitions, each having slightly
#different transition matrices 

F=120
N=4
T=12

stateNames = seq(10000,10000*N,10000)
for(k in 1:N){
	stateNames[k]=toString(stateNames[k])
}


#base transition matrix that makes it less likely to transition to states that are far away
normal = 0
basematrix = matrix(0L,nrow=N,ncol=N)
for (i in 1:N){
	vec = seq(0,0,length.out=N)
	normal = 0
	for (j in 1:N){
		p = 1/2*1/(abs(i-j)+1)
		normal = normal+p
		vec[j]=p
	}
	basematrix[i,]=vec*1/normal
}


#vary the transition matrix for each firm -- STARTING HERE
matrixList = list()

#Randomly create market trend
mtrend = runif(1)

for(k in 1:F){
	matrixList[[k]]=matrix(0L,nrow=N,ncol=N)

	#Does the firm tend to do better or worse?
	#The better the market trend is, the more likely the firm is to do well
	f = runif(1)
	ftrend = runif(1)
	ftrend = sign(mtrend-f)*ftrend*1/8


	for(i in 1:N){
		vec = seq(0,0,length.out=N)
		normal = 0
		for (j in 1:N){
			
			if (i-j<0){
				#if the trend is positive the probability of getting into a higher sales class in increased
				pInc = rnorm(1,ftrend,2*abs(ftrend))*basematrix[i,j]+basematrix[i,j]
			} else if (i-j>0) {
				#if the trend is positive the probability of falling into a lower sales class is decreased
				pInc = rnorm(1,(-1)*ftrend,2*abs(ftrend))*basematrix[i,j]+basematrix[i,j]
			} else {
				pInc = rnorm(1,1/16,1/8)*basematrix[i,j]+basematrix[i,j]
			}			
			normal = normal+pInc
			vec[j]=pInc
		}
		matrixList[[k]][i,]=vec*1/normal
	}

}

#vary the transition matrix for each firm -- ENDING HERE


#Randomized starting vector assigns firms to the states (state boundaries)

n = matrix(0,nrow=T,ncol=N+1)
n[1,] = c(0,sort(round(runif(N-1)*F)),F)


#Simulate the transitions and create contingiency tables -- STARTING HERE
firms = array(rep(c(1:F),T),dim=c(F,T))
states = matrix(0,nrow=F,ncol=T)

contin = list()

for(k in 1:(T-1)){
	contin[[k]]=matrix(0,ncol=N,nrow=N)
}

statesNew = list()
#Create T-1 contingiency tables
for (k in 1:(T-1)){

#Create Matrix for new state assignments
for(i in 1:N){
statesNew[[i]]=vector()
}

for( j in 1:N){
	for (i in firms[(n[k,j]+1):n[k,j+1],k]){
		#For all firms in this state
		states[i,k]=j
		#Generate random number to decide new state
 		x1 = runif(1)
 		

		#Find the index of the first cumulative P where x1<=P using the transition matrix for this firm
		ns = 1		
		while(x1>sum(matrixList[[i]][j,1:ns])){
			ns = ns+1
		}

		states[i,k+1]=ns

		statesNew[[ns]]=c(statesNew[[ns]],i)
		contin[[k]][j,ns] = contin[[k]][j,ns]+1
	}

}

#assign new state boundaries
for(i in 2:N){
	length = 0
	for(j in 1:(i-1)){
		length = length + length(statesNew[[j]])
	}
	n[k+1,i]=length
	
}
n[k+1,N+1]=F

#assign new order of the firms
newFirms=vector()
for(i in 1:N){
	newFirms = c(newFirms,statesNew[[i]])
}
firms[,k+1] = newFirms

}

#Simulate the transitions and create contingiency tables -- ENDING HERE


#Simulate employment levels -- STARTING HERE

L = matrix(0,ncol=T,nrow=F)

L[,1]=round(1000*abs(rnorm(F,0,1)))

Av = (N+1)/2

for(k in 2:T){
	for(f in 1:F){
		s = states[f,k-1]

		#At time k-1, firm f was in state s

		if(s-Av>0){
			newemp = abs(round(rnorm(1,(s-Av)*L[f,k-1]/50,2)))
		} else if (s-Av<0){
			newemp = -abs(round(rnorm(1,(s-Av)*L[f,k-1]/50,2)))
		} else{
			newemp = round(rnorm(1,0,2))		
		}

		L[f,k]=L[f,k-1]+newemp

	}
}

#Simulate employment levels -- ENDING HERE