% implements continuous time moel 

clear all
close all

dt = 0.01; 
T = 500; % make bigger than plot for numerical integrals to infinity
tv = (0:dt:T)';
Nt = size(tv,1);
% general setup for Lucas curves


% single exponential 

sig = 0.2; 
kap = 0.2; 
a = 1; 
etai = 0.3; 
i0 = 0.01; 
% p0 given pi0
pi0 = -0.005; 
p0 = 1/a*(pi0+sig*kap*(pi0-i0));
disp('assumption, p0 given pi0 in percent')
disp(100*p0)
% pi0 given p0. Better in the end, revert if you redo everythng. 
%p0 = -0.005; 
%pi0 = (a*p0 + sig*kap*i0)/(1+sig*kap);

s0 = 1; 
r = 0.02; % used for fiscal calculations. 
phi = 1.1; 


tv = (0:dt:T)';
Nt = size(tv,1);
pit = zeros(Nt,1);
pt = zeros(Nt,1);
xt = zeros(Nt,1);

ipterm = zeros(Nt,1);
ipiterm = zeros(Nt,1);

it = i0*exp(-etai*tv);

% Dynamic IS 

it = i0*exp(-etai*tv);

pit = (1+sig*kap).*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2)*(pi0+sig*kap*(pi0-i0))...
+sig*kap*a*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2).*(1-exp(-etai*tv))/etai*i0...
+sig*kap./(sig*kap+exp(-a*tv)).*it;
pt = (1+sig*kap)./(sig*kap+exp(-a*tv))*p0 + sig*kap./(sig*kap+exp(-a*tv)).*(1-exp(-etai.*tv))/etai*i0;
xt = exp(-a*tv).*pt/kap;

ut = it - phi*pit; 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*xt,'-k','LineWidth',2)
plot(tv,100*ut,'-m','LineWidth',2)
plot(tv,0*tv,'-k')
axis([0 10 -inf inf])
xlabel('Time','FontSize',16)
ylabel('Percent','FontSize',16)
legend('i','\pi','x','u','fontsize',16,'location','best')
print -depsc2 irdynamicct.eps

% s from p0 
% new method using i - pi 

int1 = sum(exp(-(a+r)*tv)*(1+sig*kap)./((sig*kap+exp(-a*tv)).^2)*dt );
if ~(etai==0);
    int2 = sum(exp(-(a+r)*tv)*sig*kap*a.*(1-exp(-etai*tv))./((sig*kap+exp(-a*tv)).^2*etai)*dt);
elseif (etai==0);
    int2 = sum(exp(-(a+r)*tv)*sig*kap*a./((sig*kap+exp(-a*tv)).^2)*dt);
end;
int3 = sum(exp(-(a+r+etai)*tv)./(sig*kap+exp(-a*tv))*dt);
sr = -p0*(1+a*int1)-i0*int2+i0*int3; 

% new formula from long term, should be equivalent
int4 = sum(sig*kap*exp(-(r+etai)*tv)./(sig*kap+exp(-a*tv))*dt);
sr2 = -p0*(1+a*int1)-i0*int2-i0*int4 + i0*(1/(r+etai)); 

sr1 = sum(exp(-r.*tv).*(it-pit)*dt)-p0;

disp('r*100 s0/r, new s0/r direct ');
disp([r*100 sr sr2 sr1 ]);

omeg = [0 0.0026 0.1 0.3 1 100];
sr4 = -p0*(1+a*int1)-i0*int2-i0*int4 + i0*(1/(r+etai)) - i0./(etai+r+omeg); 
disp('long term debt omega and surplus ')
disp(omeg)
disp(sr4)



%%%%%%%%%% 
% New double exponential model


sig = 0.2; 
kap = 0.2; 
a = 1; 
b = 2;
etai = 0.3; 
i0 = 0.01; 
r = 0.02; % used for fiscal calculations. 
phi = 1.1; 
%p0 = [-0.005 ]; % inherit p0 from last one 

tv = (0:dt:T)';
Nt = size(tv,1);
pit = zeros(Nt,1);
pt = zeros(Nt,1);
xt = zeros(Nt,1);
it = i0*exp(-etai*tv);

% plot a 

jv = (0:0.01:5)';
alph = (exp(-a*jv)-exp(-b*jv))/(1/a-1/b);

jv2 = (2:0.01:5)';
alph2 = (exp(-a*jv2)-exp(-b*jv2))/(1/a-1/b);

figure;

hold on
area(jv2,alph2,'facecolor','g');
plot(jv,alph,'-r','linewidth',2)
plot(jv,a*exp(-a*jv),'-b','linewidth',2)
xlabel('Time lag j','FontSize',16)
ylabel('Weight \alpha_j','FontSize',16)
text(0.5,0.7,'Exponential','FontSize',16,'Color','b')
text(1.2,0.5,'Double exponential','FontSize',16,'Color','r')
text(2.5,0.2,'\Phi(t) = \int_{j=t}^\infty \alpha_j dj','FontSize',16)
%legend('Double exponential','Exponential','fontsize',16)

print -depsc2 doublexp.eps

% calculate 

term1 = ((b*exp(-a*tv)-a*exp(-b*tv))/(b-a)+sig*kap);
term2 = a*b*(exp(-a*tv)-exp(-b*tv))/(b-a);

pit = (1+sig*kap)*term2./(term1.^2)*p0 + ...
    sig*kap*term2./(term1.^2).*(1-etai.^tv)/etai*i0 + ...
    sig*kap./term1.*it; 

pt = (1+sig*kap)./term1*p0 + sig*kap./term1.*(1-exp(-etai*tv))/etai*i0;
ptlong=[zeros(1/dt,size(pt,2));pt];
piat = ptlong(1/dt+1:end,:)-ptlong(1:end-1/dt,:);

ut = it - phi*pit; 
xt = (b*exp(-a*tv)-a*exp(-b*tv))/(b-a).*pt/kap;

%sig2 = 0.4; attempt to show price puzzle, not worth it. 
%kap2 = 0.4 
%term1a = ((b*exp(-a*tv)-a*exp(-b*tv))/(b-a)+sig2*kap2);

%pit2 = (1+sig2*kap2)*term2./(term1a.^2)*p0 + ...
%    sig2*kap2*term2./(term1a.^2).*(1-etai.^tv)/etai*i0 + ...
%    sig2*kap2./term1a.*it; 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*xt(:,1),'-k','LineWidth',2)
plot(tv,100*pit(:,1),'-r','LineWidth',2)
plot(tv,100*piat(:,1),'-m','LineWidth',2)
plot(tv(1:1/dt:end),100*piat(1:1/dt:end,1),'vm','MarkerFaceColor','m')

%plot(tv,100*pit2(:,1),'--r','LineWidth',2)
plot(tv,0*tv,'-k')
axis([0 10 -inf inf])
xlabel('Time','FontSize',16)
ylabel('Percent','FontSize',16)
legend('i','x','\pi','p_t-p_{t-1}','fontsize',16,'location','best')
print -depsc2 doublexpresp.eps

%
% display multiple equilibria -- different pi0

pi0 = -[ 0 0.25 0.5 1]/100;
p0 = 1/a*(pi0+sig*kap*(pi0-i0));% p0 given pi0
pi0 = ones(Nt,1)*pi0;
p0 = ones(Nt,1)*p0;

tv = tv*ones(1,size(pi0,2));
it = it*ones(1,size(pi0,2));

pit = (1+sig*kap).*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2).*(pi0+sig*kap*(pi0-i0))...
+sig*kap*a*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2).*(1-exp(-etai*tv))/etai*i0...
+sig*kap./(sig*kap+exp(-a*tv)).*it;
pt = (1+sig*kap)./(sig*kap+exp(-a*tv)).*p0 + sig*kap./(sig*kap+exp(-a*tv)).*(1-exp(-etai.*tv))/etai*i0;
xt = exp(-a*tv).*pt/kap;

figure;
hold on
plot(tv(:,1),100*it(:,1),'-b','LineWidth',2)
plot(tv(:,1),100*pit(:,1),'--r','LineWidth',2)
plot(tv(:,1),100*pit(:,2),'--r','LineWidth',2)
plot(tv(:,1),100*pit(:,3),'-r','LineWidth',2)
plot(tv(:,1),100*pit(:,4),'--r','LineWidth',2)
plot(tv(:,1),0*tv(:,1),'-k')
axis([0 10 -inf inf])
xlabel('Time','FontSize',16)
ylabel('Percent','FontSize',16)
legend('i',['\pi, \pi_0=' num2str(100*pi0(1,1))],['\pi, \pi_0=' num2str(100*pi0(1,2))],['\pi, \pi_0=' num2str(100*pi0(1,3))],['\pi, \pi_0=' num2str(100*pi0(1,4))],'fontsize',16,'location','best')
print -depsc2 irmulteqpi0.eps

% s from p0

int1 = sum(exp(-(a+r)*tv)*(1+sig*kap)./((sig*kap+exp(-a*tv)).^2)*dt );
if ~(etai==0);
    int2 = sum(exp(-(a+r)*tv)*sig*kap*a.*(1-exp(-etai*tv))./((sig*kap+exp(-a*tv)).^2*etai)*dt);
elseif (etai==0);
    int2 = sum(exp(-(a+r)*tv)*sig*kap*a./((sig*kap+exp(-a*tv)).^2)*dt);
end;
int3 = sum(exp(-(a+r+etai)*tv)./(sig*kap+exp(-a*tv))*dt);

sr = -p0(1,:).*(1+a*int1)-i0*int2+i0*int3; 
sr1 = sum(exp(-r.*tv).*(it-pit)*dt)-p0(1,:);
        
  
disp('surplus to hit initial inflation assumptions. pi0, integral, direct')
disp(pi0(1,:)*100)
disp(sr) 
disp(sr1)

%
% display multiple equilibria --interest rates 
%

pi0 = -0.005;
i0 =   [ 0  0.01 0.05 0.01];
etai = [0.3  0.3 0.3    0 ]; % don't use very small eta i. (1-e^-etai t)/etai isn't happy
tv = tv(:,1);
it = exp(-tv*etai).*i0 - [ 0 0 0.03 0];
i0 = it(1,:);
p0 = 1/a*(pi0+sig*kap*(pi0-i0));% p0 given pi

for indx = 1:size(i0,2);
    if etai(indx) > 0;
        pit(:,indx) = (1+sig*kap).*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2)*(pi0+sig*kap*(pi0-i0(indx)))...
            +sig*kap*a*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2).*(1-exp(-tv(:,1)*etai(indx)))./etai(indx).*i0(indx)...
            +sig*kap./(sig*kap+exp(-a*tv)).*it(:,indx);
        pt(:,indx) = (1+sig*kap)./(sig*kap+exp(-a*tv))*p0(indx) ...
            + sig*kap./(sig*kap+exp(-a*tv)).*(1-exp(-tv*etai(indx)))./etai(indx).*i0(indx);
        xt(:,indx) = exp(-a*tv).*pt(:,indx)/kap;
    elseif etai(indx)==0;
        pit(:,indx) = (1+sig*kap).*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2)*(pi0+sig*kap*(pi0-i0(indx)))...
            +sig*kap*a*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2).*i0(indx)...
            +sig*kap./(sig*kap+exp(-a*tv)).*it(:,indx);
        pt(:,indx) = (1+sig*kap)./(sig*kap+exp(-a*tv))*p0(indx) ...
            + sig*kap./(sig*kap+exp(-a*tv)).*i0(indx);
        xt(:,indx) = exp(-a*tv).*pt(:,indx)/kap;
    end
end  

figure;
hold on
plot(tv(:,1),100*it(:,1),'--b','LineWidth',2)
plot(tv(:,1),100*it(:,2),'-b','LineWidth',2)
plot(tv(:,1),100*it(:,3),'--b','LineWidth',2)
plot(tv(:,1),100*it(:,4),'--b','LineWidth',2)

plot(tv(:,1),100*pit(:,1),'--r','LineWidth',2)
plot(tv(:,1),100*pit(:,2),'-r','LineWidth',2)
plot(tv(:,1),100*pit(:,3),'--r','LineWidth',2)
plot(tv(:,1),100*pit(:,4),'--r','LineWidth',2)
plot(tv(:,1),0*tv(:,1),'-k')
axis([0 10 -inf inf])
xlabel('Time','FontSize',16)
ylabel('Percent','FontSize',16)
%legend('i','\pi, \pi_0=0','\pi, \pi_0=-0.5%','\pi, \pi_0=-1%','\pi, \pi_0=-1.5%','fontsize',16,'location','best')
print -depsc2 irmulteqi0.eps

% s from p0 

int1 = sum(exp(-(a+r)*tv)*(1+sig*kap)./((sig*kap+exp(-a*tv)).^2)*dt );
int2 = ones(1,size(etai,2));
int3 = ones(1,size(etai,2));
for indx = 1:size(int2,2);
    if ~(etai(indx)==0);
        int2(indx) = sum(exp(-(a+r)*tv)*sig*kap*a.*(1-exp(-etai(indx)*tv))./((sig*kap+exp(-a*tv)).^2*etai(indx))*dt);
    elseif (etai(indx)==0);
        int2(indx) = sum(exp(-(a+r)*tv)*sig*kap*a./((sig*kap+exp(-a*tv)).^2)*dt);
    end;
    int3(indx) = sum(exp(-(a+r+etai(indx))*tv)./(sig*kap+exp(-a*tv))*dt);
end;

sr = -p0(1,:).*(1+a*int1)-i0.*int2+i0.*int3; 
sr1 = sum(exp(-r.*tv).*(it-pit)*dt)-p0;
        
disp('surplus for different interest rate paths.i eta pi0 s/r and direct')
disp('do not use integral for permanent negative shock. Integral assumes AR(1)')
disp(it(1,:)*100)
disp(etai);
disp(pi0(1,:)*100)
disp(sr) 
disp(sr1)


% %%%%%%%%%%
% Calvo pricing
 
% classic Taylor Rule
 
etai = 0.3; 
phi = 1.1; 
r = 0.02; 
kap = etai/sqrt(2); 
sig = etai/sqrt(2); % hard wired to produce pi = -i
i0 = 0.01;
u0 = (1+sig*kap*phi/(etai*(r+etai)-sig*kap))*i0 ; 

disp('Calvo parameters etai phi r kap sig i0 u0')
disp([etai phi r kap sig i0 u0])

 %    x    pi         u
A = [ 0 sig*(phi-1) sig ;
     -kap  r           0 ; 
     0     0          -etai];

[Q L] = eig(A); 
Q1 = inv(Q);
Lv = diag(L); 
indpls = find(real(Lv) > 0);
Q1pls = Q1(indpls,:);
Qpls = Q(:,indpls);
Lvpls = Lv(indpls);
indmin = find(real(Lv) <= 0);
Q1min = Q1(indmin,:);
Qmin = Q(:,indmin);
Lvmin = Lv(indmin);

indknown0 = 3; % index of known x0 
xknown0 = [u0]; 
indfree0 = [1 2]; %index of unknown x0
xu = -inv(Q1pls(:,indfree0))*Q1pls(:,indknown0)*xknown0;
x0 = zeros(size(A,1),1);
x0(indknown0) = xknown0;
x0(indfree0) = xu; 
zmin0 = Q1min*x0;
zmint = exp(tv*Lvmin').*(ones(size(tv,1),1)*zmin0');
xt = zmint*Qmin';

yt = xt(:,1);
pit = xt(:,2);
ut = xt(:,3);
it = phi*pit + ut; 

% analytic

lam1 = (r + sqrt(r^2-4*sig*kap*(phi-1)))/2;
lam2 = (r - sqrt(r^2-4*sig*kap*(phi-1)))/2;
ut2 = u0*exp(-etai*tv);
pit2 = sig*kap/(lam1-lam2)*(1/(lam1+etai)-1/(lam2+etai))*ut2;
pit3 = -sig*kap/((lam1+etai)*(lam2+etai))*ut2;
pit4 = -sig*kap/(sig*kap*(phi-1)+etai*r + etai^2)*ut2;

% figure 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
%plot(tv,100*pit2,'-r','LineWidth',2)
%plot(tv,100*pit3,'-r','LineWidth',2)
%plot(tv,100*pit4,'-r','LineWidth',2)% checks, they're all =

plot(tv,100*yt,'-k','LineWidth',2)
plot(tv,100*ut,'-m','LineWidth',2)
%plot(tv,ut2,'mv','LineWidth',2)
plot(tv,0*tv,'-k')
legend('i','\pi','x','u','fontsize',16)
text(1,0.9,'i','color','b','fontsize',18)
text(1,-0.5,'\pi','color','r','fontsize',18)
text(2,1.4,'u','color','m','fontsize',18)
text(2,-0.8,'x','color','k','fontsize',18)

axis([0 10 -inf inf])
xlabel('Time','fontsize',16)
ylabel('Percent','fontsize',18)
print -depsc2 taylor1.eps
print -dpng taylor1.png

disp('standard calvo')
disp('s0/r for peramanent surplus taylor rule. Computed and direct')
disp(-etai/(sig*kap)*pit(1));
sr1 = sum(exp(-r.*tv).*(it-pit)*dt);
disp(sr1)


% alternative solution with pi0 = 0. 

T = 100 ; 
etai=0.3;
rho = 0.02; 
pi0 = 0;
i0 = 0.01; 

tv = (0:dt:T)';
Nt = size(tv,1);
pit = zeros(Nt,1);
xt = zeros(Nt,1);
Pt = zeros(Nt,1);
it = i0*exp(-etai*tv);
pit(1) = pi0;

lam1 = (-rho + sqrt(rho^2+4*sig*kap))/2;
lam2 = (rho + sqrt(rho^2+4*sig*kap))/2;

% analytical solution
pi0 = 0;
pit = pi0 * exp(-lam1*tv)... 
+ i0 * lam1*lam2/((etai-lam1)*(lam2+etai)) * (exp(-lam1*tv)-exp(-etai*tv));
xt = 1/kap*( pi0*(rho*lam1)*(i0*exp(-lam1*tv)))...
    + i0*lam1*lam2/((etai-lam1)*(etai+lam2)) *...
    ((rho+lam1)*exp(-lam1*tv)-(rho+etai)*exp(-etai*tv));

phi = 1.1; 
vt = it - phi*pit; 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*xt,'-k','LineWidth',2)
plot(tv,100*vt,'-m','linewidth',2)
plot(tv,0*tv,'-k')
legend('i','\pi','x','u','fontsize',16)
xlabel('Time','fontsize',16)
ylabel('Percent','FontSize',16)
axis([0 10 -inf inf]);
print -depsc2 taylor2.eps
print -dpng taylor2.png

disp('calvo with pi0=0')
disp('s0/r for peramanent surplus taylor rule.  direct')
sr1 = sum(exp(-r.*tv).*(it-pit)*dt);
disp(sr1)

pit2 = pit;
xt2 = xt; 
vt2 = vt; % save 


% Calvo, find zero surplus solution 

pi0 = i0/((etai-lam1)*(etai+lam2))*...
    ((r+lam1)*etai*(etai+lam2-lam1)/(r+etai)-lam1*lam2);

pit = pi0 * exp(-lam1*tv)... 
+ i0 * lam1*lam2/((etai-lam1)*(lam2+etai)) * (exp(-lam1*tv)-exp(-etai*tv));
xt = 1/kap*( pi0*(rho*lam1)*(i0*exp(-lam1*tv)))...
    + i0*lam1*lam2/((etai-lam1)*(etai+lam2)) *...
    ((rho+lam1)*exp(-lam1*tv)-(rho+etai)*exp(-etai*tv));

phi = 1.1; 
vt = it - phi*pit; 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*xt,'-k','LineWidth',2)
plot(tv,100*vt,'-m','linewidth',2)
plot(tv,0*tv,'-k')
legend('i','\pi','x','u','fontsize',16)
xlabel('Time','fontsize',16)
ylabel('Percent','FontSize',16)
axis([0 10 -inf inf]);
text(2,0.7,'i','color','b','fontsize',18)
text(1,0.5,'\pi','color','r','fontsize',18)
text(1.5,0.25,'u','color','m','fontsize',18)
text(5,0.07,'x','color','k','fontsize',18)
print -depsc2 taylor3.eps
print -dpng taylor3.png

disp('calvo with s=0')
disp('s0/r for peramanent surplus taylor rule.  direct')
sr1 = sum(exp(-r.*tv).*(it-pit)*dt);
disp(sr1)

% both on one graph? 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit2,'--r','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*vt2,'--m','linewidth',2)
plot(tv,100*vt,'-m','linewidth',2)

plot(tv,0*tv,'-k')
legend('i','\pi, \pi_0=0', '\pi, s=0', 'u, \pi_0=0','u, s=0','fontsize',16)
xlabel('Time','fontsize',16)
ylabel('Percent','FontSize',16)
axis([0 10 -inf inf]);
%text(2,0.7,'i','color','b','fontsize',18)
%text(1,0.5,'\pi','color','r','fontsize',18)
%text(1.5,0.25,'u','color','m','fontsize',18)
%text(5,0.07,'x','color','k','fontsize',18)
print -depsc2 taylor4.eps


% Mixed IS 


sig = 0.2; 
kap = 0.2; 
a = 1; 
etai = 0.3; 
i0 = 0.01; 
pi0 = -0.005; 
s0 = 1; 
r = 0.02; % used for fiscal calculations. 

m = [ 1 0.5 ];
Nm = size(m,2);
pit = zeros(Nt,Nm);
xt = zeros(Nt,Nm);
ipterm = zeros(Nt,Nm);
ipiterm = zeros(Nt,Nm);
it = i0*exp(-etai*tv);

mm = ones(Nt,1)*m; 
tm = tv*ones(1,Nm);

% interest rate term 
for ind = (1:Nt);
    t = tv(ind);
    if ind==1; 
        ipterm(ind,:)  = 0;
        ipiterm(ind,:) = 0; 
    else
        ipterm(ind,:)  = sig*kap./((m*exp(-a*t)+sig*kap).^((1-m)/(a*m)+1)).* ...
        sum((exp(-a*tv(1:ind-1))+sig*kap).^((1-m)./(a*m)).*(it(1:ind-1)*ones(1,Nm))*dt);   
        ipiterm(ind,:)  = sum( ((mm(1:ind-1,:).*exp(-a*tm(1:ind-1,:))+sig*kap)./...
              (mm(1:ind-1,:).*exp(-a*t)+sig*kap)).^...
                 ((1-mm(1:ind-1,:))./(a*mm(1:ind-1,:)))...
                 .*(it(1:ind-1)*ones(1,Nm))*dt);  
    end;
end;


pit2 = exp(-a*tm).*((mm+sig*kap)./(mm.*exp(-a*tm)+sig*kap)).^((1-mm)./(a*mm)+2).*...
    ((mm*pi0+sig*kap*(pi0-i0))./(mm+sig*kap))+...
    (1-mm*(1-a)).*((sig*kap*exp(-a*tm))./(mm.*exp(-a*tm)+sig*kap).^2).*ipiterm+...
    (sig*kap)./(mm.*exp(-a.*tm)+sig*kap).*it; 
pt2 = ((mm+sig*kap)./(mm.*exp(-a*tm)+sig*kap)).^((1-mm*(1-a))./(a*mm)).*...
    ((mm*pi0 + sig*kap*(pi0-i0))./(1-mm*(1-a)))+ipterm;
xt2 = exp(-a*tm).*pt2/kap;

if 0
    figure;
    hold on
    plot(tv,it,'-b','LineWidth',2)
    plot(tv,pt2,'-r','LineWidth',2)
    plot(tv,xt2,'-k','LineWidth',2)
    plot(tv,0*tv,'-k')
    legend('i','p','x')
end

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit2(:,1),'--r','LineWidth',2)
plot(tv,100*pit2(:,2),'-r','LineWidth',2)
plot(tv,100*xt2(:,1),'--k','LineWidth',2)
plot(tv,100*xt2(:,2),'-k','LineWidth',2)
plot(tv,0*tv,'-k')
legend('i',['\pi, m = ' num2str(m(1))],['\pi, m =' num2str(m(2))],...
    ['x, m = ' num2str(m(1))], ['x, m = ' num2str(m(2))],'fontsize',16,'location','best')
ylabel('Percent','fontsize',16)
xlabel('Time','fontsize',16)
axis([0 10 -inf inf])
print -depsc2 mixed.eps

disp('downweighted IS ')
disp('s0 for peramanent surplus taylor rule.  direct')
sr1 = sum(exp(-r.*tv).*(it-pit2)*dt);
disp(sr1)


% AR(1) habits 

th = 0.7; 
b = 1; 

pi0 = -0.005; 
tv = (0:dt:T)';
Nt = size(tv,1);
pit = zeros(Nt,1);
xt = zeros(Nt,1);
Pt = zeros(Nt,1);
pt = zeros(Nt,1);
it = i0*exp(-etai*tv);
p0 = 1/(a+th*b)*(pi0+sig*kap*(pi0-i0));
pt(1) = p0; 
pit(1) = pi0;

% solve via differential equation 
for ind = (1:Nt-1);
    t = tv(ind);
    dpdt = ((a+th*b)*exp(-a*t)*pt(ind)...
        -th*b*Pt(ind) + sig*kap*it(ind))/(exp(-a*t)+sig*kap);
    pt(ind+1) = pt(ind) + dpdt*dt;
    Pt(ind+1) = Pt(ind) + b*(exp(-a*t)*pt(ind)-Pt(ind))*dt; % P here is kappa X in the paper
    if t>0;
        pit(ind) = dpdt;
    end;
end;
xt = exp(-a*tv).*pt/kap; 

sr1 = sum(exp(-r.*tv).*(it-pit)*dt)-p0;
disp('habits surplus')
disp(sr1)


if 0
figure;
hold on
plot(tv,it,'-b','LineWidth',2)
plot(tv,pt,'-r','LineWidth',2)
plot(tv,xt,'-k','LineWidth',2)
plot(tv,Pt,'-m','LineWidth',2)
plot(tv,0*tv,'-k')
legend('i','p','x','P')
end

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*xt,'-k','LineWidth',2)
plot(tv,0*tv,'-k')
legend('i','\pi','x','fontsize',16,'location','best')
xlabel('Time','fontsize',16);
ylabel('Percent','fontsize',16);
axis([0 10 -inf inf])
print -depsc2 habit.eps

% AR(1) habits Calvo Pricing 
 

th = 0.7; 
b = 1; 
pi0 = -0.005; 

etai = 0.3; 
phi = 1.1; 
r = 0.02; 
kap = 0.2; 
sig = 0.2; % hard wired to produce pi = -i
i0 = 0.01;
 %    X    x    pi     i
A = [    -b     b        0     0 ;
      -th*b  th*b     -sig   sig ;   
          0  -kap        r     0 ; 
          0     0        0 -etai];

[Q L] = eig(A); 
Q1 = inv(Q);
Lv = diag(L); 
indpls = find(real(Lv) > 0);
Q1pls = Q1(indpls,:);
Qpls = Q(:,indpls);
Lvpls = Lv(indpls);
indmin = find(real(Lv) <= 0);
Q1min = Q1(indmin,:);
Qmin = Q(:,indmin);
Lvmin = Lv(indmin);

indknown0 = [1 3 4]; % index of known x0, inflation  
xknown0 = [0 pi0 i0]'; 
indfree0 =  [2]; %index of unknown x0
xu = -inv(Q1pls(:,indfree0))*Q1pls(:,indknown0)*xknown0;
x0 = zeros(size(A,1),1);
x0(indknown0) = xknown0;
x0(indfree0) = xu; 
zmin0 = Q1min*x0;
zmint = exp(tv*Lvmin').*(ones(size(tv,1),1)*zmin0');
xt = zmint*Qmin';

Xt = xt(:,1);
yt = xt(:,2);
pit = xt(:,3);
it = xt(:,4);


% figure 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*yt,'-k','LineWidth',2)
plot(tv,0*tv,'-k')
legend('i','\pi','x','fontsize',16)
axis([0 10 -inf inf])
xlabel('Time','fontsize',16)
ylabel('Percent','fontsize',16)
print -depsc2 habitcalvo.eps

sr1 = sum(exp(-r.*tv).*(it-pit)*dt);
disp('habits + calvo surplus')
disp(sr1)

%%%%%%%%%%%%%
% Long term debt


sig = 0.2; 
kap = 0.2; 
a = 1; 
etai = 0.05; 
i0 = 0.01; 
pi0 = -0.005; 
s0 = 1; 
r = 0.02; % used for fiscal calculations. 
phi = 1.1; 


tv = (0:dt:T)';
Nt = size(tv,1);
pit = zeros(Nt,1);
pt = zeros(Nt,1);
xt = zeros(Nt,1);

ipterm = zeros(Nt,1);
ipiterm = zeros(Nt,1);

it = i0*exp(-etai*tv);

% Dynamic IS 

it = i0*exp(-etai*tv);
p0 = 1/a*(pi0+sig*kap*(pi0-i0));% p0 given pi0

pit = (1+sig*kap).*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2)*(pi0+sig*kap*(pi0-i0))...
+sig*kap*a*exp(-a*tv)./((sig*kap+exp(-a*tv)).^2).*(1-exp(-etai*tv))/etai*i0...
+sig*kap./(sig*kap+exp(-a*tv)).*it;
pt = (1+sig*kap)./(sig*kap+exp(-a*tv))*p0 + sig*kap./(sig*kap+exp(-a*tv)).*(1-exp(-etai.*tv))/etai*i0;
xt = exp(-a*tv).*pt/kap;

ut = it - phi*pit; 

figure;
hold on
plot(tv,100*it,'-b','LineWidth',2)
plot(tv,100*pit,'-r','LineWidth',2)
plot(tv,100*xt,'-k','LineWidth',2)
plot(tv,100*ut,'-m','LineWidth',2)
plot(tv,0*tv,'-k')
axis([0 10 -inf inf])
xlabel('Time','FontSize',16)
ylabel('Percent','FontSize',16)
legend('i','\pi','x','u','fontsize',16,'location','best')
print -depsc2 longdebtct.eps

% s from p0 
% new method using i - pi 

int1 = sum(exp(-(a+r)*tv)*(1+sig*kap)./((sig*kap+exp(-a*tv)).^2)*dt );
if ~(etai==0);
    int2 = sum(exp(-(a+r)*tv)*sig*kap*a.*(1-exp(-etai*tv))./((sig*kap+exp(-a*tv)).^2*etai)*dt);
elseif (etai==0);
    int2 = sum(exp(-(a+r)*tv)*sig*kap*a./((sig*kap+exp(-a*tv)).^2)*dt);
end;
int3 = sum(exp(-(a+r+etai)*tv)./(sig*kap+exp(-a*tv))*dt);
sr = -p0*(1+a*int1)-i0*int2+i0*int3; 

% new formula from long term, should be equivalent
int4 = sum(sig*kap*exp(-(r+etai)*tv)./(sig*kap+exp(-a*tv))*dt);
sr2 = -p0*(1+a*int1)-i0*int2-i0*int4 + i0*(1/(r+etai)); 

sr1 = sum(exp(-r.*tv).*(it-pit)*dt)-p0;

disp('r*100 s0/r, new s0/r direct ');
disp([r*100 sr sr2 sr1 ]);

omegzero = i0/sr2-r-etai; 

omeg = [0 omegzero 0.1 0.3 1 100];
sr4 = -p0*(1+a*int1)-i0*int2-i0*int4 + i0*(1/(r+etai)) - i0./(etai+r+omeg); 
disp('long term debt')
disp(omeg)
disp(sr4)

