clc
clear all

global sigma psi chi alpha gamma beta delta kappa T gs epsilon theta phi lambda StatVars rBar yBar IfSave

Vars=[];
MaxFunEvals=50000;
MaxIter=200;
options=optimset('Display','iter','MaxFunEvals',MaxFunEvals,'MaxIter',MaxIter,'Algorithm','levenberg-marquardt','TolFun',.00000000000000000000001,'TolX',.0000000000000000001);
params=[
        0.6,0.5;
        0.7,0.5;
        0.8,0.5;
        0.89,0.5
        0.9,0.5;
        0.95,0.5;
        0.89,-1.0
        ];
sz=size(params);
Ncase=sz(1);
Nparam=sz(2);
T=80;
for icase=1:Ncase
    theta=params(icase,1);
    phi=params(icase,2);
    Vars=[Vars;params(icase,:),T,zeros(1,15-Nparam)];
    sigma=.4;
    psi=1.54;
    chi=0.334;
    gamma=1.103;
    beta=0.95^0.25;
    delta=0.1/4;
    kappa=8;
    alpha=0.7;
    pi=.7^0.25;
    epsilon=5;
   
    lambda=0;
    
    %  Static solution with g=0
    
    x(1:3,1)=[1;1;3];
    res=NKStatModel(x);
    MaxFunEvals=50000;
    MaxIter=200;
    options=optimset('Display','iter','MaxFunEvals',MaxFunEvals,'MaxIter',MaxIter,'Algorithm','levenberg-marquardt','TolFun',.00000000000000000000001,'TolX',.0000000000000000001);
    xNew=fsolve(@NKStatModel,x,options);
    resNew=NKStatModel(xNew);
    c=xNew(1);
    h=xNew(2);
    K=xNew(3);
    g=0;
    k=K-c;
    muc=c^(-1/sigma)*(1-chi*(1-1/sigma)*h^(1+1/psi));
    y=h^alpha*k^(1-alpha);
    yBar=y;
    cBar=c;
    mu=epsilon/(epsilon-1);
    muh=(-chi*(1+1/psi)*c^(1-1/sigma)-gamma)*h^(1/psi);
    w=alpha*h^(alpha-1)*k^(1-alpha)/mu;
    r=(1-alpha)*h^alpha*k^(-alpha)/mu-delta;
    rBar=r-phi-lambda;
    b=r+delta;
    StatVars=[y,1,h,r,w,b,1,c,k,0]; % [y,p,h,r,w,b,q,c,k,g]
    Vars=[Vars;StatVars,zeros(1,3),1/mu,0,0];  %  [y,p,h,r,w,b,q,c,k,g,ms,ptildes,Ds,zs,As,Bs]
    
    %  Dynamic solution
    
    %  The unknowns are values of y, p, h, w, b, q, and c for t=2:T,
    %  r for t=3:T, and
    %  k for t=2:T-1
    %  Boundary conditions are for p and k in t=1 and k in T
    
    yLoc=1;
    pLoc=yLoc+T-1;
    hLoc=pLoc+T-1;
    rLoc=hLoc+T-1;
    wLoc=rLoc+T-2;
    bLoc=wLoc+T-1;
    qLoc=bLoc+T-1;
    cLoc=qLoc+T-1;
    kLoc=cLoc+T-1;
    x(yLoc:yLoc+T-2,1)=y;
    x(pLoc:pLoc+T-2)=1;  %   initial value for p is 1
    x(hLoc:hLoc+T-2)=h;
    x(rLoc:rLoc+T-3)=r;
    x(wLoc:wLoc+T-2)=w;
    x(bLoc:bLoc+T-2)=b;
    x(qLoc:qLoc+T-2)=1;  %  initial value for Tobin's q
    x(cLoc:cLoc+T-2)=c;
    x(kLoc:kLoc+T-3)=k;
    g=.01;
    for t=2:T
        gs(t,1)=g;
        g=g*pi;
    end
    IfSave=1;
    disc=NKACCModel(x);
    MaxFunEvals=50000;
    MaxIter=30;
    options=optimset('Display','iter','MaxFunEvals',MaxFunEvals,'MaxIter',MaxIter,'Algorithm','levenberg-marquardt','TolFun',.00000000000000000000001,'TolX',.0000000000000000001);
    IfSave=0;
    xN=fsolve(@NKACCModel,x,options);
    IfSave=1;
    Check=NKACCModel(xN);
    
    %  Extract variables
    
    load NKACCModel AllVars
    Vars=[Vars;AllVars];
    mults(icase,1)=(AllVars(2,1)-yBar)/gs(2); % y multiplier
    mults(icase,2)=(AllVars(2,8)-cBar)/gs(2);  % c multiplier
    mults(icase,3)=log(AllVars(2,14)*mu)/log(AllVars(2,1)/yBar);  %  omega
end

mults=[mults,params];
xlswrite('Results.xlsx',mults,'NKmults');
Params=[sigma,psi,chi,gamma,beta,delta,kappa,alpha,pi,epsilon,phi,lambda]';
xlswrite('Results.xlsx',Params,'NKParams');
% xlswrite('Results.xlsx',Vars,'NKVars');
