%This code computes the wage dynamics of the Stock-Flow paper for several
%values of k.
clear all
close all
format long

%Set parameter values
x=1;                            % match productivity 
b=0.4;                          % opportunity cost of employment
lambda=0.45;                    % Probability of an UE flow (Shimer, 2005b) : rate at which vacancies enter the job centre
delta=0.02;                    % Probability of a Job to Job Transition (Nagypal, 2008) : rate at which vacancies leave the job centre
mu=0.02;                        % Probability of an EU flow (Shimer, 2005a) : rate at which workers experience a re-set shock
r=0.0041;                       % time discount rate.
T=0.000001:0.05:30.05;          % time since last visit
K=[0,1,2,5,10,20];    % number of firms left in the job centre since the last visit
phi=zeros(1,length(T));
phip=zeros(1,length(T));
Prob=zeros(length(T),length(K));
W=zeros(length(T),length(K));
Prob2=zeros(length(T),length(K));

%WAGE EQUATION 
C=(x-b)/(r+mu+delta+lambda);
for k=1:length(K)
    for t=1:length(T)
    phi(t)=(lambda/delta)*(1-exp(-delta*T(t)));
    phip(t)=lambda*exp(-delta*T(t));
    Prob(t,k)=(exp(-phi(t)))*(1-exp(-delta*T(t)))^K(k);
    W(t,k)=x-Prob(t,k)*((r+mu+delta)-((phip(t)/phi(t))*(K(k)-phi(t))))*C;
    Prob2(t,k)=Prob(t,k)*(phip(t)/(phi(t)^2))*(phip(t)*((K(k)-phi(t))^2)-K(k)*(delta*phi(t)+phip(t))+delta*(phi(t)^2));
    end
end

%Graphs

plot(T,W(:,1),T,W(:,2))
axis([0,max(T),(b-0.05),1.3])
pause
plot(T,W(:,1),T,W(:,2),T,W(:,3))
axis([0,max(T),(b-0.05),1.3])
pause
plot(T,W(:,1),T,W(:,2),T,W(:,3),T,W(:,4))
axis([0,max(T),(b-0.05),1.3])
pause
plot(T,W(:,1),T,W(:,2),T,W(:,3),T,W(:,4),T,W(:,5))
axis([0,max(T),(b-0.05),1.3])
pause
plot(T,W(:,1),T,W(:,2),T,W(:,3),T,W(:,4),T,W(:,5),T,W(:,6))
axis([0,max(T),(b-0.05),1.3])