clc;
clear;
close all;

%global beta delta alpha cplusKp Kbarfix pr Vold K count;
beta = 0.99;
delta = 0.025;
alpha = 0.36;
tolerance=.01;
maxits=300;
lowbnd=0.01;
upbnd = 25;

pr = [0.8507    0.1229    0.5833    0.0938
    0.1159    0.8361    0.0313    0.3500
    0.0243    0.0021    0.2917    0.0313
    0.0091    0.0389    0.0938    0.5250]';

nK=300;             %number of indicies in our state space
K = zeros(1,nK)';  %populating the state space so that number of possible 
for i = 1:nK+1;      %workers is a function of (1-delta).  This way we ensure 
    K(i) = 25*(1-delta)^(nK+1 - i); %that we are in the state space in the following period
end;
nKbar=10;
Kbar = (1:1:nKbar);

%Kprime(i,j)= K(I);   

Vold=ones(nKbar,nK+1,4);                      %Starting values of value function
%Notice the updated value function will only have ONE maximand for each
%state combination. That's why the second dimension is just 1.
Vnew=ones(nKbar,nK+1,4);                      %Will update after each iteration
%EVold=0;                                        %order V(Kbar,K,z,e)
%Vtemp=zeros(nKbar,nK,4);                     %will popolate with values for each choice 

dif=1;
its=0;
%loop for a given good state, employed
ltilda=1;
lbar=.99*1;  %in good state so unemp is 1%
e = [1 0];
z = [1.01 .99];
while dif>tolerance && its<maxits
    for i=1:nKbar;         %loop through average capital
        count = 0; 
        for k = 1:1:length(e)
         for j = 1:1:length(z) 
            count = count + 1;
            lbar=.99-.03*(z(j)==.99); %check this
            w = (1-alpha)*z(j)*(Kbar(i)/lbar)^alpha;
            r = alpha*z(j)*(Kbar(i)/lbar)^(alpha-1); 
             for l=1:nK+1;                     %for each possible individual capital holdings
                if (z(j)== 1.01)
                    Kbarprime = exp(.095+.962*log(Kbar(i)));
                end
                if (z(j)== 0.99)
                    Kbarprime = exp(.085+.965*log(Kbar(i)));
                end
                Kbarfix=round(Kbarprime);
                
                cplusKp = r*K(l) + w*lbar*e(k) + (1-delta)*K(l);
                % Feed in values.  Then, Find the kbarprime that maximizes the
                % function.
                %f = @(Kprime)valuefn(Kprime,cplusKp,Kbarfix,pr,Vold);
                %interpolate the value function off grid
                ZI = interp2(X,Y,Z,XI,YI)
                [KprimeMax,V] = fminbnd(@(Kprime) valuefn(Kprime,K,count,beta,cplusKp,Kbarfix,pr,Vold),lowbnd,upbnd);
                
                % Vtemp = -V since we find the minumum
                Vnew(i,l,count) = -V;
                
             end;
            %The line above calculates the maximum among possible values of
            %k'
            %Kprime(i,l)= K(I);          %optimal number of people in firm in next period for each i,j
            %Vtemp=zeros(nK,1);          %Reset Vtemp for next iteration
         end;   
        end;
    end;
    Vnewinterp = interp2(Vnew);
    Vold=Vnewinterp;
    %dif=max(max(abs(Vnew-Vold)));
    %Vold=Vnew;
    %its=its+1;
    display('Aggy loves life.')
    display(Vnew)
    display(KprimeMax);
end;

if its==maxits
    display('VFI: Max Iterations reached');
end;


display (Vold')
display (Kprime')
Vold(:,1) = Vold(:,2);
Kprime(:,1) = Kprime(:,2);
set(0,'DefaultAxesColorOrder',[0 0 0],'DefaultAxesLineStyleOrder','-|-.|--|:')

figure()
plot(K(2:end,:), Vold')
title('Graph 2a: Value function for low persistence, no firing ')
legend('low price', 'medium price','high price')
xlabel('N, Current Number of Employees')