function [eta,gdp,yrs,ytrend]=transitiondynamics_1shock(eta0,mutildeH,mutildeL,sigmaH,sigmaL,pbar,qbar,death,NumGridPoints,alpha,s_old,s_new,g_old,Tshock,yrs_data,eta_data,figname,CaseName,Tyears,dt,PL_percentile);
% transitiondynamics_1shock.m
%
%  Basic code for transition dynamics solving the KFE
%  when there is a single shock to phiH, sigmaH, tau, or pBar, for example.

disp ' ';
N=NumGridPoints; % Number of Grid Points for x, e.g. 1000
if g_old>.5; g_old=g_old/100; end;

% Create the initial SS distribution of x ==> f0
%  -- Find the pbar0 that delivers eta0
mm= mutildeH / sigmaH^2;
pbar0=1/2*sigmaH^2*( (1/eta0+mm)^2 - mm^2 ) - death
[f0,g0,eta0_out,A,x,dx]=tdstationary(mutildeH,mutildeL,sigmaH,sigmaL,pbar0,qbar,death,N,PL_percentile);
fprintf(' eta0 (input)         = %5.3f\n',eta0);
fprintf(' eta0 (output,solved) = %5.3f\n',eta0_out);
%if abs(eta0-eta0_out)> 1e-3; disp 'eta0 is not eta0_out??? Stopping...'; keyboard; end;
mm= mutildeL / sigmaL^2;
xiL=-mm+sqrt(mm^2+2*(death)/sigmaL^2);
fprintf(' etaL (formula)            = %5.3f\n',1/xiL);
disp ' ';

X0 = sum( x.*f0*dx ); % mean of the initial distribution
yInit=1*X0^alpha*(1-s_old);
y(1)=yInit;


% Find the new SS distribution ==> fSS
[fSS,gSS,etaSS,A_new,x,dx]=tdstationary(mutildeH,mutildeL,sigmaH,sigmaL,pbar,qbar,death,N,PL_percentile);

% Check with the formula
mm= mutildeH / sigmaH^2;
xiH=-mm+sqrt(mm^2+2*(pbar+death)/sigmaH^2);
fprintf(' etaSS (formula)         = %5.3f\n',1/xiH);
fprintf(' etaSS (output, solved)  = %5.3f\n',etaSS);
mm= mutildeL / sigmaL^2;
xiL=-mm+sqrt(mm^2+2*(death)/sigmaL^2);
fprintf(' etaL (formula)            = %5.3f\n',1/xiL);


T=Tyears/dt;
[eta,f_t]=transdyn(A_new,g0,Tyears,dt,x,dx,PL_percentile);

% gdp
g_new=g_old*s_new/s_old;
G=(1+g_new*ones(T,1)*dt); % Growth factor is cumprod(G(1:(t-1))
G=[1; G]; % Initialize growth factor

X = sum( mult(f_t,x)*dx )'; % mean of the distribution
y=cumprod(G).*X.^alpha*(1-s_new);
y=[X(1)^alpha*(1-s_old); y];
gdp=y/y(1)*200;
gdp=trimr(gdp,0,2);
disp ' ';

yrs=Tshock+(1:T)'*dt;
oldyrs=(Tshock-10:Tshock)';
ytrend=gdp(1)*(1+g_old).^(yrs-(Tshock+1));
oldgdp=gdp(1)*(1+g_old).^(oldyrs-(Tshock));

yrs=[oldyrs; yrs];
gdp=[oldgdp; gdp];
ytrend=[oldgdp; ytrend];
eta=[eta0_out*ones(size(oldyrs)); eta];


definecolors;
LW=4;
figure(1); figsetup; makefigwide
if exist('yrs_data');
    plot(yrs_data,eta_data,'-','Color',myred,'LineWidth',2);
    plot(yrs_data,eta_data,'o','Color',myred,'LineWidth',2);
end;
plot(yrs,etaSS*ones(size(yrs)),'--','Color',mygreen,'LineWidth',1);
plot(yrs,eta,'-','Color',myblue,'LineWidth',LW);
%plot([1970 1975.1],[eta0_out eta0_out],'-','Color',myblue,'LineWidth',LW);
chadfig2('Year','Pareto Inequality, \eta',1,0);
print('-depsc',[figname '_' CaseName '_eta']);

%tt=(Tshock-10/dt):min(T,Tshock+60/dt);
tt=(1:length(yrs));
figure(2); figsetup; makefigwide
%plot([oldyrs; yrs],log([oldgdp; ytrend]),'--','Color',mygreen,'LineWidth',1);
plot(yrs,log(ytrend),'--','Color',mygreen,'LineWidth',1);
plot(yrs,log(gdp),'-','Color',myblue,'LineWidth',LW);
%plot(oldyrs,log(oldgdp),'-','Color',myblue,'LineWidth',LW);
vstep=100;
v1=floor(min(gdp(tt))/vstep)*vstep;
v2=ceil(max(gdp(tt))/vstep)*vstep;
vals=vstep*2.^(0:ceil(log2(v2/vstep))); %v1:vstep:v2;
relabelaxis(log(vals),num2str(vals'),'y');
chadfig2('Year','GDP per person',1,0);
print('-depsc',[figname '_' CaseName '_gdp']);


disp ' '; disp ' ';
if exist('ParamToShock');
    disp(['Results from Shocking ' ParamToShock]);
    ParamToShockValues
end;


data=[yrs eta gdp ytrend gdp./ytrend];
tt=[1:(1/dt):(30/dt) (35/dt):(5/dt):T]';
tle='Year Eta GDP Trend Ratio';
cshow(' ',data(tt,:),'%8.0f %8.4f %8.0f %8.0f %8.3f',tle);





function [f,g,eta,A,x,dx]=tdstationary(mutildeH,mutildeL,sigmaH,sigmaL,pbar,qbar,death,N,PL_percentile);

% tdstationary -- computes the "A" matrix for KFE
%    as well as the stationary distribution given
%    mutildeH/L, sigmaH/L, pbar, and death 
%
%   g := stack of [gH gL]
%   f := the density gH+gL


%i_entry = round(N/4);
xmin = -2;
xmax = 15; %8;
x = linspace(xmin,xmax,N)';
[x0 i_entry] = min(abs(x));
dx = (xmax-xmin)/(N-1);
dx2 = dx^2;
sH = mutildeH*ones(N,1);
sL = mutildeL*ones(N,1);

% Construct the "A" transition matrix
% A is the transition matrix for the KFE, which Moll calls A* or AT (transpose)
%  ==> df = A*f

% This is Chad's version
XH =  sigmaH^2/(2*dx2)*ones(N,1);  % The n-2 coefficients
YH =  sH/dx - sigmaH^2/dx2;        % The n-1 coefficients
ZH =  -sH/dx + sigmaH^2/(2*dx2);   % The  n  coefficients
XL =  sigmaL^2/(2*dx2)*ones(N,1);  % The n-2 coefficients
YL =  sL/dx - sigmaL^2/dx2;        % The n-1 coefficients
ZL =  -sL/dx + sigmaL^2/(2*dx2);   % The  n  coefficients


AH=spdiags(YH,0,N,N)+spdiags(XH(2:N),-1,N,N)+spdiags([0;ZH(1:N-1)],1,N,N);
AH(N,N) = YH(N) + ZH(N);
AH(1,1) = XH(1)+ YH(1);
AL=spdiags(YL,0,N,N)+spdiags(XL(2:N),-1,N,N)+spdiags([0;ZL(1:N-1)],1,N,N);
AL(N,N) = YL(N) + ZL(N);
AL(1,1) = XL(1)+ YL(1);

% Switching states and death
Aswitch = [-(death+pbar)*speye(N) sparse(N,N); pbar*speye(N) -death*speye(N)];

% Entry: qbar -> H, 1-qbar -> L -- coming from both gH and gL
Aentry=sparse(2*N,2*N);
Aentry(i_entry,:)=qbar*death;
Aentry(N+i_entry,:)=(1-qbar)*death;

% Putting it all together
%A0=[AH sparse(N,N); sparse(N,N) AL]'; % Ben
A0=[AH sparse(N,N); sparse(N,N) AL];   % Chad
A = A0 + Aswitch + Aentry;

% Solve 0=A*f to get the stationary distribution
%  Note: need to fix one value, otherwise matrix is singular
N2=2*N;
b = zeros(N2,1);
i_fix = 1;
b(i_fix)=.02; %.1;
row = [zeros(1,i_fix-1),1,zeros(1,N2-i_fix)];
AA = A;
AA(i_fix,:) = row;

% Solve for the stationary distribution
% g := stack of [gH gL]
% f := the density gH+gL
g = AA\b;
g_sum = g'*ones(N2,1)*dx;
g = g./g_sum;
f=g(1:N)+g((N+1):end); % Add the gH and gL densities together


% Compute Pareto exponent in the tail
%p = 0.99; % Top 1 percent
p = PL_percentile; % e.g. Top 0.1 percent
F=cumsum(f*dx);
[val index] = min(abs(F - p));
PL=-delchange(log(f))/dx;
PL1=PL(index);
eta=1/PL1

% % Use S(.1)/S(1) to estimate eta
% % Recall f is the density for log(income)
% F=cumsum(f*dx);
% [val i01] = min(abs(F - 0.99));
% [val i001] = min(abs(F - 0.999));
% S01 =sum(exp(x(i01:end)).*f(i01:end)*dx);
% S001 =sum(exp(x(i001:end)).*f(i001:end)*dx);
% eta=1+log10(S001/S01);


% % Checking the Pareto index
% % Note well: as NumGridPoints gets large (e.g. 15000), the index converges
% % to the stationary value at the top. For coarser grids, there is some
% % discrepancy. Also depends on how far up you look...
% i1=find(x>6);
% i2=find(x>13);
% tt=(i1:i2);
% figure(1); clf;
% plot(x(i1:i2),1./PL(i1:i2));

% figure(2); clf;
% tt=index+(-50:50);
% plot(x(tt),1./PL(tt));

% disp 'Stopping in tdstationary';
% keyboard;



function [eta,f]=transdyn(A,g0,Tyears,dt,x,dx,PL_percentile);

% Transition dynamics
%   Iterates the discretized KFE to solve for the time
%   path of the distribution
T=Tyears/dt;
N2=size(A,1);
N=N2/2;

g=zeros(N2,T+1);
g(:,1) = g0;
for t=1:T;
    g(:,t+1) = (speye(N2,N2) - dt*A)\g(:,t);
end

ff=g(1:N,:)+g((N+1):end,:); % Add the gH and gL densities together
                            
f_sum=sum(ff)*dx;
f = div(ff,f_sum);

p = PL_percentile; % e.g. 0.999 for Top 0.1 percent
F=cumsum(f*dx);
for t=1:T
    [val index] = min(abs(F(:,t) - p));
    PL1_t(t) = -(log(f(index+1,t))-log(f(index,t)))/dx;
end
eta=1./PL1_t';

% % Use S(.1)/S(1) to estimate eta
% % Recall f is the density for log(income)
% F=cumsum(f*dx);
% for t=1:T;
%     [val i01] = min(abs(F(:,t) - 0.99));
%     [val i001] = min(abs(F(:,t) - 0.999));
%     S01 =sum(exp(x(i01:end)).*f(i01:end,t)*dx);
%     S001 =sum(exp(x(i001:end)).*f(i001:end,t)*dx);
%     eta(t)=1+log10(S001/S01);
% end;
% eta=eta';

%keyboard
