%%%%% Problem set 2

clear all

%%%%% QUESTION 1: STOCHASTIC DIFFERENCE EQUATIONS

display('OUTPUT FOR QUESTION #1')

x0    = 0;
sigma = 1;

Ts     = [100,500];  NT = 2;
As     = [0.5,0.99]; NA = 2;
Bs     = [0,1];      NB = 2;

X_short = [];
X_long  = [];
means   = [];
stddevs = [];


for tt = 1:2,
    
    T = Ts(tt);

    Z = randn(T,1);
    X = zeros(T,1);
    
    for aa=1:2,
        
    a = As(aa);    
        
    for bb=1:2,    
    
    b = Bs(bb);    
        
    for t=1:T-1,       
        X(t+1) = a*X(t)+b+sigma*Z(t+1);
    end    
    
        if T==100,
            X_short = [X_short,X];
            means   = [means,b/(1-a)];
            stddevs = [stddevs,sigma/sqrt(1-(a^2))];
        else
            X_long  = [X_long,X];
            means   = [means,b/(1-a)];
            stddevs = [stddevs,sigma/sqrt(1-(a^2))];
        end
        
end
end
end

time100 = cumsum(ones(100,1));
time500 = cumsum(ones(500,1));

figure(1)
subplot(4,2,1)
plot(time100,X_short(:,1),time100,means(1)*ones(100,1),'k--',...
    time100,(means(1)-2*stddevs(1))*ones(100,1),'r--',time100,(means(1)+2*stddevs(1))*ones(100,1),'r--')
ylabel(['a = 0.50, b = 0.0, T = 100'])
subplot(4,2,2)
plot(time100,X_short(:,2),time100,means(2)*ones(100,1),'k--',...
    time100,(means(2)-2*stddevs(2))*ones(100,1),'r--',time100,(means(2)+2*stddevs(2))*ones(100,1),'r--')
ylabel(['a = 0.50, b = 1.0, T = 100'])
subplot(4,2,3)
plot(time100,X_short(:,3),time100,means(1)*ones(100,1),'k--',...
    time100,(means(3)-2*stddevs(3))*ones(100,1),'r--',time100,(means(3)+2*stddevs(3))*ones(100,1),'r--')
ylabel(['a = 0.99, b = 0.0, T = 100'])
subplot(4,2,4)
plot(time100,X_short(:,4),time100,means(4)*ones(100,1),'k--',...
    time100,(means(4)-2*stddevs(4))*ones(100,1),'r--',time100,(means(4)+2*stddevs(4))*ones(100,1),'r--')
ylabel(['a = 0.99, b = 1.0, T = 100'])

subplot(4,2,5)
plot(time500,X_long(:,1),time500,means(5)*ones(500,1),'k--',...
    time500,(means(5)-2*stddevs(5))*ones(500,1),'r--',time500,(means(5)+2*stddevs(5))*ones(500,1),'r--')
ylabel(['a = 0.50, b = 0.0, T = 500'])
subplot(4,2,6)
plot(time500,X_long(:,2),time500,means(6)*ones(500,1),'k--',...
    time500,(means(6)-2*stddevs(6))*ones(500,1),'r--',time500,(means(6)+2*stddevs(6))*ones(500,1),'r--')
ylabel(['a = 0.50, b = 1.0, T = 500'])
subplot(4,2,7)
plot(time500,X_long(:,3),time500,means(7)*ones(500,1),'k--',...
    time500,(means(7)-2*stddevs(7))*ones(500,1),'r--',time500,(means(7)+2*stddevs(7))*ones(500,1),'r--')
ylabel(['a = 0.99, b = 0.0, T = 500'])
subplot(4,2,8)
plot(time500,X_long(:,4),time500,means(8)*ones(500,1),'k--',...
    time500,(means(8)-2*stddevs(8))*ones(500,1),'r--',time500,(means(8)+2*stddevs(8))*ones(500,1),'r--')
ylabel(['a = 0.99, b = 1.0, T = 500'])

%break %if you want to stop here!! 

clear all

%%%%% QUESTION 2: MARKOV CHAINS

display('OUTPUT FOR QUESTION #2')

mu    = 0.02;
sigma = 0.04;

x     = [mu+sigma;mu-sigma];

pi0   = [1,0];


pp    = [0.1,0.5,0.9];
TT    = [50,100,1000];

means   = [];
stddevs = [];
autos   = [];
ps      = [];
Ts      = [];

for i = 1:3,
    
p     = pp(i);

P     = [p,1-p;1-p,p];

for j = 1:3,

T     = TT(j);

[X,states] = simulate_markov(x,P,pi0,T);

xbar = mean(X);
s2   = std(X);

%%%%%  To do the autocorrelation, we need to lag X
%%%%%  We lose one degree of freedom

lagX = X(1:T-1);
X    = X(2:T);

rho  = corrcoef(X,lagX); %% this produces a 2-by-2 matrix of correlation coefficients, we want the off-diagonal element
rho  = rho(1,2);

means   = [means,xbar];
stddevs = [stddevs,s2];
autos   = [autos,rho];
ps      = [ps,p];
Ts      = [Ts,T];

end
end

ps
Ts
means
stddevs
autos

%break %if you want to stop here!! 

clear all

%%%%% QUESTION 3: PREFERENCES OVER MARKOV CHAINS

display('OUTPUT FOR QUESTION #3')

c   = [1;5];
pi0 = [0.5,0.5];


beta  = 0.95;
gamma = 2.5;
P1    = eye(2,2);
P2    = 0.5*ones(2,2);

u   = (1/(1-gamma))*c.^(1-gamma);

v1    = inv(eye(2,2)-beta*P1)*u;
v2    = inv(eye(2,2)-beta*P2)*u;

display('gamma = 2.5')

V1    = pi0*v1
V2    = pi0*v2

gamma = 4.0;
P1    = eye(2,2);
P2    = 0.5*ones(2,2);

u   = (1/(1-gamma))*c.^(1-gamma);

v1    = inv(eye(2,2)-beta*P1)*u;
v2    = inv(eye(2,2)-beta*P2)*u;

display('gamma = 4.0')

V1    = pi0*v1
V2    = pi0*v2

%break %if you want to stop here!! 

clear all

%%%%% QUESTION 4: COUPLED MARKOV CHAINS

display('OUTPUT FOR QUESTION #4')

%%%%% Aggregate state

mu    = 0.02;
sigma = 0.04;
g     = [mu+sigma;mu-sigma];
Q     = [0.99,0.01;0.01,0.99];
pig   = [1,0];

%%%%% Individual employment state

e     = [1;0];
Ph    = [0.99,0.01;0.99,0.01];
Pl    = [0.50,0.50;0.50,0.50];
pie   = [1,0];

%%%%% Coupled chain

x     = [g(1),e(1);g(1),e(2);g(2),e(1);g(2),e(2)];

P     = [Ph*Q(1,1),Pl*Q(1,2);Ph*Q(2,1),Pl*Q(2,2)];

pi0   = [1,0,0,0];

%%%%% Stationary distribution of Px chain

[V,D] = eig(P');
d     = diag(D);
[smallest,index] = min(abs(d-1)); %% find a unit eigenvalue
v     = V(:,index);

pibar = (v/sum(v))'

%%%%% Simulate chain

%%%%% Simulation (4 subplots)

T = 250;

for j = 1:4,

[chain,states] = simulate_markov(x,P,pi0,T);

G = chain(1,:);
E = chain(2,:);

W = zeros(size(E));

for i = 1:length(E),
    if E(i)==1,
        W(i)=0.10;
    else
        W(i)=0.08;
    end
end

ts = cumsum(ones(T,1));

figure(2)
subplot(2,2,j)
plot(ts,G,'bo-',ts,W,'rx-')
axis([1 T -0.04 0.14])

end