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