proc iml;

/* generate a realization of a markov chain, given the transition matrix and initial state */

n = 1200;			/* total length of chain */
n0 = 100;			/* number to trim from front */
nt = 100;			/*    and back of chain */

s = 1;			/* initial state */

/* transition matrix, ij element = P[move to state i from state j] */
b = {0.7 0.1 0.1 0.1, 0.1 0.7 0.1 0.1, 0.1 0.1 0.7 0.1, 0.1 0.1 0.1 0.7};

k = 5;		/* number of stages in the transition matrix */

/* all elements of the transition matrix entered in row order, each row */
/*    corresponds to a different environment */

alla = {0 0 0 0 0               0    0.422 0.065 0.174 0.071     0 0.244 0.455 0.174 0.286     0 0.022 0.351 0.628 0.286     0 0 0 0.047 0.571,
        0 0 0.001 0.013 0.041   0.50 0.808 0.166 0.095 0.103     0 0.239 0.621 0.175 0.138     0 0.047 0.187 0.574 0.207     0 0 0.006 0.052 0.586,
		0 0 0.001 0.011 0.033   0.455 0.481 0.167 0.100 0.206    0 0.176 0.583 0.345 0.294     0 0.030 0.119 0.437 0.559     0 0 0.005 0.026 0.353,
		0 0 0.013 0.286 0.869   0.750 0.588 0.101 0.107 0.048    0 0.320 0.592 0.266 0.048     0 0.039 0.223 0.469 0.286     0 0 0.018 0.170 0.571};

/* print out each transition matrix just to make sure they were entered correctly */

A0 = J(k,k,0);
do i = 1 to nrow(alla);
  a = shape(alla[i,],k,k);
  print a;
  a0 = a0+a;
  end;

amean = a0/nrow(alla);
call eigen(eval, evec, amean);
maxeval = loc(eval[,1]=max(eval[,1]));		/* locate the maximum real eigenvalue */	
lmean = eval[maxeval,1];
loglmean = log(lmean);

print lmean;
print loglmean;

call eigen(eval, evec, b);					/* calculate eigenvect and val of B */
maxeval = loc(eval[,1]=max(eval[,1]));		/* locate the maximum real eigenvalue */	
u = evec[,maxeval]/sum(evec[,maxeval]);		/* stable state distribution */
print b;
print eval;
print u;

cb = j(nrow(b),ncol(b),0);			/* need cum sum within columns to generate states */
do i = 1 to ncol(b);
  cb[,i] = cusum(b[,i]);
  end;

c = j(n,1,0);		/* space to store the sequence of environments */

do i = 1 to n;	/* generate each state sequentially */
  s = 1+nrow(b)-sum(uniform(0) < cb[,s]);	
     /* s'th (current state) column of b (and cb) describe distribution of new states */
	 /*   generate a discrete r.v. with the given probabilities */
  c[i] = s;		/* then store new state in the chain */
  end;

bit = t(c[1:20]);
print bit;

u0 = j(k,1,1)/k;		/* initial stage distribution is equally divided among categories */
u = j(k, n, 0);			/* space to store all stage vectors */
u[,1] = u0;				/*   put initial one in place */
l = j(n,1,0);			/* and create space to store the lambda */

do i = 2 to n;
  at =shape(alla[c[i],],k,k);	/* extract At for the environment this year */
  au = at*u[,i-1];			/* do this multiplication only once */

  l[i] = sum(au);			/* store this lambda */
  u[,i] = au/l[i];			/* and stable size distribution */
  end;

goodl = log(l[(n0+1):(n-nt)]);
meanlogl = sum(goodl)/(n-nt-n0);
selogl = sqrt(ssq(goodl-meanlogl)/(n-nt-n0-1));

print meanlogl;
print selogl;

quit;