%DISCRETE DISTRIBUTIONS -- MARKOV CHAINS
%Stick wiuth foodlevel example
%food level is 'low', 'medium' or 'high'
%f(1), f(2), f(3)-- a vector of ENVIRONMENTAL STATES

%
%Model 2. Now we suppose that there is a CORRELATION between
%food levels from one year to the next -- but only the next
%this is called a Markov process -- tomorrow depends directly on today
%but not on yesterday



%How to make a model? The idea is to look at TRANSITIONS between
%states -- we ask, if state is 1 (low food) this year, what are the
%probabilities of staying in 1 (low food again), or jumping up to
%state 2 (med food) or even to state 3 (hi food).
%For Example. If food depends on
%weather, env quality, etc, it may be reasonable to suppose that
%next year will be similar to this year unless thee is a large
%change. In this case the most likely thing will be to go from
%state 1 to state 1 etc

%IN GENERAL we must specify TRANSITION PROBABILITIES
%c(i,j) = Prob{state at time t+1 is j GIVEN that state at time t is i}

%**************note that the convention in stochastic process models
%is the opposite of what we are used to in demography!!!!*************
% so in this program rows and columns are opposite to what we do
%but when you make your matrix later today you shoud follow the demographic
%convention that the colum state is now(t)and the row state is later (t+1)
%***************************************************************************


%In the present case there are 3 states, so c ill be a 3 by 3 matrix
%Because we must go to SOME state next year starting from ANY state
%this year, we have

%Sum Over j { c(i,j) } = 1

%

%EXAMPLE -- Assuming that tomorrow is more like today than something new

%TRY

c = [0.6 0.2 0.2;0.2 0.6 0.2; 0.2 0.2 0.6]

%PROBABILITY VECTORS -- suppose that at time t=1 the probability of
%state 1 is q0(1), etc. Then the prob of being in State 1 at t=2 is
% q1(2) = q0(1)*c(1,1) + q0(2)*c(2,1) + q0(3)*c(3,1)

%THEREFORE

%if prob distr of states at time 0 is q0 = (q0(1), q0(2), q0(3))
%then prob distr of states at time 0 is q1 = (q1(1), q1(2), q1(3))
% and  q1 = q0*c;  -- The q's MUST be ROW vectors -- you need to be
%careful about rows & columns in Matlab

%EXAMPLE

q0=[0.5 0.2 0.3]; % a row

q1 = q0*c %at time 1

q3 = q0*c*c*c  %at time 3

%
%You will note that at time t (for any integer t)
% qt = q0* c^t

%

%So try t = 10 & t = 11, then t = 12
%GO AHEAD

%YOU will see that q11 & q12 are ALMOST THE SAME
%IN OTHER WORDS -- THE MARKOV PROCESS CONVERGES TO
%AN EQUILIBRIUM PROBABILILTY DISTRIBUTION OF STATES

%  WE WILL CALL THIS qstar
%  OBVIOUSLY qstar = qstar*c --

%  THIS mans that qstar is a LEFT eigenvector of c with eigenvalue 1

% Code for finding qstar---why??
[vec,ev] = eig(c');
[evmax,loc] =max(abs(diag(ev)));

qstar = vec(:,loc);
qstar = qstar/sum(qstar);


%SIMULATION FROM A Markov chain

%This is really not hard -- remember that in a Markov chain
%at each time t there is a DISCRETE set of states
%so if we know the prob distr at time t, it is just a discrete
%distribution and we can just do what we did in EX3 --

%THUS we will simulate a STATE at t=1
%NEXT, one particular row of matrix c tells us the distribution
%of states at t=2 GIVEN the state we have just found

%And so on, ad infinitum
%So to SIMULATE the STATES for

T=100;

EnvState = zeros(1,T); %we will fill this up with states
%start with the q0 we had earlier and find
t=1;
probvec = q0;
psum = cumsum(probvec);
X = rand(1,1);
j = min(find(psum > X));
EnvState(1) = j;

%So at time t=1 we are in the state j
%which means that at t=2 we will be moving from state j


t=t+1;
%the following line was corrected relative to a previous version
probvec = c(j,:);% find the row j...
                 % previous version had c(:,j), which finds the column
                 % it worked in the test c, because it was symmetric
psum = cumsum(probvec);
X = rand(1,1);
j = min(find(psum > X));
EnvState(1) = j;

%And so on

%Easier to write a LOOP function
T=100;
EnvState = zeros(1,T);
probvec = q0;
X = rand(1,T); %faster to get all random numbers at one go

for time = 1:T
   psum = cumsum(probvec);
j = min(find(psum > X(time)));
EnvState(time) = j;
%the following line is corrected relative to previous version
probvec = c(j,:);% find the row j...
                 % previous version had c(:,j), which finds the column
                 % it worked in the test c, because it was symmetric
end

EnvState

%EXPLORE


figure(1)
plot(1:T, EnvState) %Note the RUNS in the environmental state
%
%and SHARP transitions

figure(2)
hist(EnvState)


%CORRELATIONS -- I have noted that we assumed that tomorrow is like
%today in our EXAMPLE of c (construct one where tomorrow is UNLIKE today!)
%to see this it helps to look at correlations between different'
%times t
% so make

v = [EnvState(1:80); EnvState(2:81);EnvState(3:82);
EnvState(4:83);EnvState(5:84);EnvState(6:85);
EnvState(7:86);EnvState(8:87);EnvState(9:88);
EnvState(10:89)];
v = v';

%Here first col is observations starting at time 1
%second col is obs from time 2 , etc
%so we have LAGGED the observations
%now

EnCor = corrcoef(v); %read help corrcoef --

%

%EnCor is a correlatiuon structure over time

% so EnCor(1,5) is the correlation between Environmental
%States that are 4 = (5-1) years apart


figure(3)
plot(1:10, EnCor(1,:))



%Note that the Correlation goes to ZERO as time separation
%Increases BUT there appears to be
%some correlation at some years --
%Oscillations in the correlation are related
%to the other eigenvalues of c -- look at eig(c)
%all eigenvalues are real
%which maens that there will be NO oscillations
%the plot shows a spurious oscillation that will disappear
%if we use more samples