nstart=0.0065;
nend=-0.0121;
pargammastart=0.8180;
pargammaend=0.9601;

pardelta=0.1;
parbeta=0.9673;
parsigma=0.5;
Rast=1.04;
gast=0.006553;
alphar=0.3;
alphaw=0.5;
c=0.167;
rhoL=0.61;
omega=0.9778;
tauB=0;
zetaT=0.33;
zetaN=0.33;

options = optimset('Display', 'on');
%% find starting period condition
options = optimset('Display', 'off');
Psistart=(1-omega)/(1+nstart-pargammastart);
RKstart=Rast*(1+tauB)-(1-pardelta);
DAstart=1+gast;
PiNstart=DAstart^(1-zetaN);
Pirstart=PiNstart^(1-alphar);
Piwstart=PiNstart^(1-alphaw);


[guessresult,faval,exitflag]=fsolve(@(guess)AAstequation(guess,DAstart,rhoL,RKstart,omega,Rast,alphaw,alphar,nstart,tauB,zetaT,zetaN,pardelta,pargammastart,Pirstart,Piwstart,parsigma,parbeta),[Rast*(1+0.5*tauB),0.5],options);
if exitflag <1
    %indicate the SS computation was not sucessful; this would also be detected by Dynare
    %setting the indicator here shows how to use this functionality to
    %filter out parameter draws
    check=1; %set failure indicator
    return; %return without updating steady states
end
Rstart=guessresult(1);
TNwTLTstart=guessresult(2);
xirstart=1-pargammastart*parbeta^(parsigma)*(Rstart/Pirstart)^(parsigma-1);
findxiw=@(xiw)(1/xiw-1-parbeta^(parsigma)*(Rstart/Piwstart*(omega+(1-omega)*(xirstart/xiw)^(1/(1-parsigma))))^(parsigma-1)/xiw);
[xiwstart,faval,exitflag]=fzero(findxiw,xirstart,options);
if exitflag <1
    %indicate the SS computation was not sucessful; this would also be detected by Dynare
    %setting the indicator here shows how to use this functionality to
    %filter out parameter draws
    check=1; %set failure indicator
    return; %return without updating steady states
end
Omegastart=omega+(1-omega)*(xirstart/xiwstart)^(1/(1-parsigma));
lambdastart=(1-omega)/(1-omega*(1-xirstart)*Rstart/DAstart/(1+nstart));
AAststart=1-(DAstart-1)/(c*TNwTLTstart);
TANwTKTstart=(zetaT*TNwTLTstart^(1-zetaT)/RKstart)^(1/(1-zetaT))*DAstart*(1+nstart);
TATwTstart=RKstart*TANwTKTstart/(DAstart*(1+nstart))/TNwTLTstart*(1-zetaT)/zetaT;
TNwTLNstart=(1-TNwTLTstart^((rhoL+1)/rhoL))^(rhoL/(rhoL+1));
TATwNstart=(TNwTLNstart/TNwTLTstart)^(1/rhoL)*TATwTstart;
TANwTKNstart=TATwNstart*TNwTLNstart*zetaN/(1-zetaN)*DAstart*(1+nstart)/RKstart;
TNwTYNstart=(TANwTKNstart/(DAstart*(1+nstart)))^(zetaN)*TNwTLNstart^(1-zetaN);
TANwTYTstart=(TANwTKTstart/(DAstart*(1+nstart)))^(zetaT)*TNwTLTstart^(1-zetaT);
TATPNstart=TATwNstart*TNwTLNstart/(1-zetaN)/TNwTYNstart;
TANwTHstart= (TATwTstart*TNwTLTstart+TATwNstart*TNwTLNstart)/(1-omega*DAstart/(Omegastart*Rstart));
TANwTBstart=(TATPNstart*TNwTYNstart-(1-alphaw)*xiwstart*TANwTHstart)/((1-alphaw)*xiwstart*(1-lambdastart)*Rstart/(DAstart*(1+nstart))+(1-alphar)*xirstart*lambdastart*Rstart/(DAstart*(1+nstart)));
TANwTCTstart=alphaw*xiwstart*((1-lambdastart)*Rstart*TANwTBstart/(DAstart*(1+nstart))+TANwTHstart)+alphar*xirstart*lambdastart*Rstart*TANwTBstart/(DAstart*(1+nstart));
TANwTBFstart=TANwTBstart-(TANwTKTstart+TANwTKNstart);
TNwTCNstart=TNwTYNstart;

%% find ending period condition
Psiend=(1-omega)/(1+nend-pargammaend);
RKend=Rast*(1+tauB)-(1-pardelta);
DAend=1+gast;
PiNend=DAend^(1-zetaN);
Pirend=PiNend^(1-alphar);
Piwend=PiNend^(1-alphaw);


[guessresult,faval,exitflag]=fsolve(@(guess)AAstequation(guess,DAend,rhoL,RKend,omega,Rast,alphaw,alphar,nend,tauB,zetaT,zetaN,pardelta,pargammaend,Pirend,Piwend,parsigma,parbeta),[Rast*(1+0.5*tauB),0.5],options);
if exitflag <1
    %indicate the SS computation was not sucessful; this would also be detected by Dynare
    %setting the indicator here shows how to use this functionality to
    %filter out parameter draws
    check=1; %set failure indicator
    return; %return without updating steady states
end
Rend=guessresult(1);
TNwTLTend=guessresult(2);
xirend=1-pargammaend*parbeta^(parsigma)*(Rend/Pirend)^(parsigma-1);
findxiw=@(xiw)(1/xiw-1-parbeta^(parsigma)*(Rend/Piwend*(omega+(1-omega)*(xirend/xiw)^(1/(1-parsigma))))^(parsigma-1)/xiw);
[xiwend,faval,exitflag]=fzero(findxiw,xirend,options);
if exitflag <1
    %indicate the SS computation was not sucessful; this would also be detected by Dynare
    %setting the indicator here shows how to use this functionality to
    %filter out parameter draws
    check=1; %set failure indicator
    return; %return without updating steady states
end
Omegaend=omega+(1-omega)*(xirend/xiwend)^(1/(1-parsigma));
lambdaend=(1-omega)/(1-omega*(1-xirend)*Rend/DAend/(1+nend));
AAstend=1-(DAend-1)/(c*TNwTLTend);
TANwTKTend=(zetaT*TNwTLTend^(1-zetaT)/RKend)^(1/(1-zetaT))*DAend*(1+nend);
TATwTend=RKend*TANwTKTend/(DAend*(1+nend))/TNwTLTend*(1-zetaT)/zetaT;
TNwTLNend=(1-TNwTLTend^((rhoL+1)/rhoL))^(rhoL/(rhoL+1));
TATwNend=(TNwTLNend/TNwTLTend)^(1/rhoL)*TATwTend;
TANwTKNend=TATwNend*TNwTLNend*zetaN/(1-zetaN)*DAend*(1+nend)/RKend;
TNwTYNend=(TANwTKNend/(DAend*(1+nend)))^(zetaN)*TNwTLNend^(1-zetaN);
TANwTYTend=(TANwTKTend/(DAend*(1+nend)))^(zetaT)*TNwTLTend^(1-zetaT);
TATPNend=TATwNend*TNwTLNend/(1-zetaN)/TNwTYNend;
TANwTHend= (TATwTend*TNwTLTend+TATwNend*TNwTLNend)/(1-omega*DAend/(Omegaend*Rend));
TANwTBend=(TATPNend*TNwTYNend-(1-alphaw)*xiwend*TANwTHend)/((1-alphaw)*xiwend*(1-lambdaend)*Rend/(DAend*(1+nend))+(1-alphar)*xirend*lambdaend*Rend/(DAend*(1+nend)));
TANwTCTend=alphaw*xiwend*((1-lambdaend)*Rend*TANwTBend/(DAend*(1+nend))+TANwTHend)+alphar*xirend*lambdaend*Rend*TANwTBend/(DAend*(1+nend));
TANwTBFend=TANwTBend-(TANwTKTend+TANwTKNend);
TNwTCNend=TNwTYNend;


%n=[repmat(nstart,1,2) nserie1to100' repmat(nend,1,901)];
%pargamma=[repmat(pargammastart,1,2) pargammaserie1to100' repmat(pargammaend,1,901)];

n=[repmat(nstart,1,3)  repmat(nend,1,1000)];
pargamma=[repmat(pargammastart,1,3)  repmat(pargammaend,1,1000)];

Psi=[repmat(Psistart,1,2) linspace(Psistart,Psiend,100) repmat(Psiend,1,901)];
RK=[repmat(RKstart,1,2) linspace(RKstart,RKend,100) repmat(RKend,1,901)];
DA=[repmat(DAstart,1,2) linspace(DAstart,DAend,100) repmat(DAend,1,901)];
PiN=[repmat(PiNstart,1,2) linspace(PiNstart,PiNend,100) repmat(PiNend,1,901)];
Pir=[repmat(Pirstart,1,2) linspace(Pirstart,Pirend,100) repmat(Pirend,1,901)];
Piw=[repmat(Piwstart,1,2) linspace(Piwstart,Piwend,100) repmat(Piwend,1,901)];
xir=[repmat(xirstart,1,2) linspace(xirstart,xirend,100) repmat(xirend,1,901)];
xiw=[repmat(xiwstart,1,2) linspace(xiwstart,xiwend,100) repmat(xiwend,1,901)];
Omega=[repmat(Omegastart,1,2) linspace(Omegastart,Omegaend,100) repmat(Omegaend,1,901)];
lambda=[repmat(lambdastart,1,2) linspace(lambdastart,lambdaend,100) repmat(lambdaend,1,901)];
TNwTLT=[repmat(TNwTLTstart,1,2) linspace(TNwTLTstart,TNwTLTend,100) repmat(TNwTLTend,1,901)];
AAst=[repmat(AAststart,1,2) linspace(AAststart,AAstend,100) repmat(AAstend,1,901)];
TANwTKT=[repmat(TANwTKTstart,1,2) linspace(TANwTKTstart,TANwTKTend,100) repmat(TANwTKTend,1,901)];
TATwT=[repmat(TATwTstart,1,2) linspace(TATwTstart,TATwTend,100) repmat(TATwTend,1,901)];
TNwTLN=[repmat(TNwTLNstart,1,2) linspace(TNwTLNstart,TNwTLNend,100) repmat(TNwTLNend,1,901)];
TATwN=[repmat(TATwNstart,1,2) linspace(TATwNstart,TATwNend,100) repmat(TATwNend,1,901)];
TANwTKN=[repmat(TANwTKNstart,1,2) linspace(TANwTKNstart,TANwTKNend,100) repmat(TANwTKNend,1,901)];
TNwTYN=[repmat(TNwTYNstart,1,2) linspace(TNwTYNstart,TNwTYNend,100) repmat(TNwTYNend,1,901)];
TANwTYT=[repmat(TANwTYTstart,1,2) linspace(TANwTYTstart,TANwTYTend,100) repmat(TANwTYTend,1,901)];
TATPN=[repmat(TATPNstart,1,2) linspace(TATPNstart,TATPNend,100) repmat(TATPNend,1,901)];
TANwTH=[repmat(TANwTHstart,1,2) linspace(TANwTHstart,TANwTHend,100) repmat(TANwTHend,1,901)];
TANwTB=[repmat(TANwTBstart,1,2) linspace(TANwTBstart,TANwTBend,100) repmat(TANwTBend,1,901)];
TANwTCT=[repmat(TANwTCTstart,1,2) linspace(TANwTCTstart,TANwTCTend,100) repmat(TANwTCTend,1,901)];
TANwTBF=[repmat(TANwTBFstart,1,2) linspace(TANwTBFstart,TANwTBFend,100) repmat(TANwTBFend,1,901)];
TNwTCN=[repmat(TNwTCNstart,1,2) linspace(TNwTCNstart,TNwTCNend,100) repmat(TNwTCNend,1,901)];
R=[repmat(Rstart,1,2) linspace(Rstart,Rend,100) repmat(Rend,1,901)];

clearvars -except n pargamma Psi RK DA PiN Pir Piw xir xiw Omega lambda TNwTLT AAst TANwTKT TATwT TNwTLN TATwN TANwTKN TANwTYT TNwTYN TATPN TANwTH TANwTB TANwTCT TANwTBF TNwTCN R;
% save initvalfile n pargamma Psi RK DA PiN Pir Piw xir xiw Omega lambda TNwTLT AAst TANwTKT TATwT TNwTLN TATwN TANwTKN TANwTYT TNwTYN TATPN TANwTH TANwTB TANwTCT TANwTBF TNwTCN;
%% function
function [result] = AAstequation(guess,DA,rhoL,RK,omega,Rast,alphaw,alphar,n,tauB,zetaT,zetaN,pardelta,pargamma,Pir,Piw,parsigma,parbeta)
    options = optimset('Display', 'off');
    result=nan(2,1);
    R=guess(1);
    TNwTLT=guess(2);
    xir=1-pargamma*parbeta^(parsigma)*(R/Pir)^(parsigma-1);

    findxiw=@(xiw)(1/xiw-1-parbeta^(parsigma)*(R/Piw*(omega+(1-omega)*(xir/xiw)^(1/(1-parsigma))))^(parsigma-1)/xiw);
    [xiw,faval,exitflag]=fzero(findxiw,xir,options);
    if exitflag <1
        %indicate the SS computation was not sucessful; this would also be detected by Dynare
        %setting the indicator here shows how to use this functionality to
        %filter out parameter draws
        check=1; %set failure indicator
        return; %return without updating steady states
    end
    Omega=omega+(1-omega)*(xir/xiw)^(1/(1-parsigma));
    lambda=(1-omega)/(1-omega*(1-xir)*R/DA/(1+n));
    TANwTKT=(zetaT*TNwTLT^(1-zetaT)/RK)^(1/(1-zetaT))*DA*(1+n);
    TATwT=RK*TANwTKT/(DA*(1+n))/TNwTLT*(1-zetaT)/zetaT;
    TNwTLN=(1-TNwTLT^((rhoL+1)/rhoL))^(rhoL/(rhoL+1));
    TATwN=(TNwTLN/TNwTLT)^(1/rhoL)*TATwT;
    TANwTKN=TATwN*TNwTLN*zetaN/(1-zetaN)*DA*(1+n)/RK;
    TNwTYN=(TANwTKN/(DA*(1+n)))^(zetaN)*TNwTLN^(1-zetaN);
    TANwTYT=(TANwTKT/(DA*(1+n)))^(zetaT)*TNwTLT^(1-zetaT);
    TATPN=TATwN*TNwTLN/(1-zetaN)/TNwTYN;
    TANwTH= (TATwT*TNwTLT+TATwN*TNwTLN)/(1-omega*DA/(Omega*R));
    TANwTB=(TATPN*TNwTYN-(1-alphaw)*xiw*TANwTH)/((1-alphaw)*xiw*(1-lambda)*R/(DA*(1+n))+(1-alphar)*xir*lambda*R/(DA*(1+n)));
    TANwTCT1=alphaw*xiw*((1-lambda)*R*TANwTB/(DA*(1+n))+TANwTH)+alphar*xir*lambda*R*TANwTB/(DA*(1+n));
    TANwTBF=TANwTB-(TANwTKT+TANwTKN);
    TANwTCT2=TANwTYT-TANwTBF+Rast*TANwTBF/(DA*(1+n))-(TANwTKN+TANwTKT)+(1-pardelta)*(TANwTKN+TANwTKT)/(DA*(1+n));
    R1=((RK+(1-pardelta))*(TANwTKN+TANwTKT)+Rast*TANwTBF)/TANwTB;
    result(1)=TANwTCT1-TANwTCT2;
    result(2)=R-R1;
end
