clc
clear all

global sigma psi chi alpha gamma beta delta kappa omega yBar T gs StatVars IfSave g

MaxFunEvals=5000;
MaxIter=10;
gCase=1;  %   1 for gradual decline of G; 2 for actual stimulus; 3 for alt stimulus
options=optimset('Display','iter','MaxFunEvals',MaxFunEvals,'MaxIter',MaxIter,'Algorithm','levenberg-marquardt','TolFun',.00000000001,'TolX',.000000001);
%  Row of params = [sigma, psi, chi, gamma, kappa, omega, gCase]
params=[
    0.4,1.54,0.334,1.103,8,0.7,1;  %  Base case
    0.4,1.54,0.334,1.103,8,0.0,1;  %  Constant markup
    0.4,1.54,0.334,1.103,0,0.7,1;  %  No adjustment cost
    0.5,1.90,0.000,1.102,8,0.7,1;  %  No complementarity
    0.4,0.50,0.334,0.617,8,0.7,1   %  Less elastic LS
    0.4,1.54,0.334,1.103,8,0.7,4   %  Knock out wealth effect
    0.4,1.54,0.334,1.103,8,1.29,2;  %  Actual stimulus
    0.4,1.54,0.334,1.103,8,1.29,3;  %  Accelerated stimulus
    0.4,1.54,0.334,1.103,8,0.7,5   %  Base case with large negative shock
    ];

sz=size(params);
Ncase=sz(1);
Nparams=sz(2);
NVars=11;
Vars=[];
Stats=[];
for icase=1:Ncase
    sigma=params(icase,1);
    psi=params(icase,2);
    chi=params(icase,3);
    gamma=params(icase,4);
    kappa=params(icase,5);
    omega=params(icase,6);
    gCase=params(icase,7);
    beta=0.95^0.25;
    delta=0.1/4;
    alpha=0.7;
    phi=.7^0.25;
    T=80;
    
    %  Static solution
    
    g=0.2;
    x=[1;1;10];
    MaxFunEvals=50000;
    MaxIter=500;
    options=optimset('Display','iter','MaxFunEvals',MaxFunEvals,'MaxIter',MaxIter,'Algorithm','levenberg-marquardt','TolFun',.00000000001,'TolX',.000000001);
    xNew=fsolve(@StatModel,x,options);
    resNew=StatModel(xNew);
    c=xNew(1);
    h=xNew(2);
    k=xNew(3);
    muc=c^(-1/sigma)*(1-chi*(1-1/sigma)*h^(1+1/psi));
    y=h^alpha*k^(1-alpha);
    yBar=y;
    cBar=c;
    muh=(-chi*(1+1/psi)*c^(1-1/sigma)-gamma)*h^(1/psi);
    w=alpha*h^(alpha-1)*k^(1-alpha);
    r=(1-alpha)*h^alpha*k^(-alpha)-delta;
    StatVars=[y,h,r,w,1,c,k,0,0]; % [y,h,r,w,q,c,k,g]
    Stats=[Stats;StatVars];
    
    %  Dynamic solution
    
    %  The unknowns are values of y, h, w, q, and c for t=2:T,
    %  k for t=2:T-1, and r for t=3:T.
    %  Boundary conditions are for k in t=1 and T
    
    yLoc=1;
    hLoc=yLoc+T-1;
    rLoc=hLoc+T-1;
    wLoc=rLoc+T-2;
    qLoc=wLoc+T-1;
    cLoc=qLoc+T-1;
    kLoc=cLoc+T-1;
    x(yLoc:yLoc+T-2,1)=y;
    x(hLoc:hLoc+T-2)=h;
    x(rLoc:rLoc+T-3)=r;
    x(wLoc:wLoc+T-2)=w;
    x(qLoc:qLoc+T-2)=1;  %  initial value for Tobin's q
    x(cLoc:cLoc+T-2)=c;
    x(kLoc:kLoc+T-3)=k;
    
    % Government purchases
    
    gsS=[0.0046;0.0046;0.0046;0.0046;0.0073;0.0073;0.0073;0.0073;0.0039;0.0039;0.0039;0.0039];
    switch gCase
        case 1  %  gradual decline
            Deltag0=0.01;
            Deltag=Deltag0;
            for t=2:T
                gs(t,1)=g+Deltag;
                Deltag=Deltag*phi;
            end
        case 2  %  actual stimulus
            gs=g+yBar*[0;gsS;zeros(T-13,1)];
        case 3 %  accelerated actual stimulus
            Deltag=1;
            lambda=.7^.25;
            for t=1:12
                gAdd(t,1)=Deltag;
                Deltag=Deltag*lambda;
            end
            gs(T,1)=g;
            gs(2:13)=g+gAdd*yBar*sum(gsS)/sum(gAdd);
        case 4  %  repayment to knock out wealth effect
            Deltag0=0.001;
            Deltag=Deltag0;
            for t=2:T
                gs(t,1)=g+Deltag;
                Deltag=Deltag*phi;
            end
            repay=beta^(-T+1)*Deltag0/(1-beta*phi);
            gs(T-1)=g-repay;
        case 5
            Deltag0=-0.1;
            Deltag=Deltag0;
            for t=2:T
                gs(t,1)=g+Deltag;
                Deltag=Deltag*phi;
            end
    end
    
    hLoc=yLoc+T-1;
    rLoc=hLoc+T-1;
    wLoc=rLoc+T-2;
    qLoc=wLoc+T-1;
    cLoc=qLoc+T-1;
    kLoc=cLoc+T-1;
    IfSave=0;
    disc=VMACCModel(x);
    MaxFunEvals=50000;
    MaxIter=20;
    options=optimset('Display','iter','MaxFunEvals',MaxFunEvals,'MaxIter',MaxIter,'Algorithm','levenberg-marquardt','TolFun',.00000000001,'TolX',.000000001);
    xN=fsolve(@VMACCModel,x,options);
    
    %  Extract variables
    
    IfSave=1;
    Check=VMACCModel(xN);
    IfSave=0;
    load VMACCModel AllVars
    Vars=[Vars;y,h,r,w,1,c,k,0,0,0,g;AllVars];
    ys=AllVars(:,1);
    cs=AllVars(:,6);
    mults(icase,1)=(ys(2)-yBar)/Deltag0;
    mults(icase,2)=(cs(2)-cBar)/Deltag0;
end
xlswrite('Results.xlsx',Vars,'VMACC');
mults=[mults,params];
xlswrite('Results.xlsx',mults,'mults');

