%zero_bound_last 
%Programs for "The New-Keynesian Liquidity Trap" JME final.  
% August 13 2017

clear all
close all

% parameters

k = 1;
sig = 1; % program is not tested for sig /=1, T/=5!
kv = [ 1 2 5 20]; %vector of values to show behavior as k changes. suffix v identifies vector versions
rho = 0.05; 
imr = 0.02; % i - r in liquidity trap
r_post = 0.02; % natural rate after all is done.
T = 5;
tauv = [ 0 0.5 1 2]' ; % vector of years to postpone interest rate rise


% derived parameters
lf = 1/2*(rho+(rho^2+4*k)^0.5); % lambdas, eigenvalues
lb = 1/2*(-rho+(rho^2+4*k)^0.5);

lfv = 1/2*(rho+(rho^2+4*kv).^0.5); % lambdas, eigenvalues; vector corresp. to kv
lbv = 1/2*(-rho+(rho^2+4*kv).^0.5);

dt = 0.01; % time interval for plots
t0 = (-5:dt:0)'; % time scale for t<0 
t1 = (0:dt:T)'; % time scale for simulation t<T
t2 = (T:dt:T+5)'; % t>T

onet0 = ones(size(t0,1),1); 
onet1 = ones(size(t1,1),1); 
onet2 = ones(size(t2,1),1); %vectors of ones for multpilications

% --------------------------------------
% x and pi for backward-stable (C=0), before and after T
% others are obtained by adding more
% or less C to it
% naming: me = C=0. w=standard (werning) 0=pi_0 = 0
% 0 1 and 2 are the three time periods
% v is a vector conformable to the vector of k for price sticky limits
% --------------------------------------

pi0me = k*sig*imr/(lf+lb) * (  (exp(-lf*(0-t0))-exp(-lf*(T-t0)))/lf);
x0me =  sig*imr/(lf+lb) * (  - lb*(exp(-lf*(0-t0))-exp(-lf*(T-t0)))/lf);

pi1me = k*sig*imr/(lf+lb) * ( (1-exp(-lb*t1))/lb + (1-exp(-lf*(T-t1)))/lf);
x1me =  sig*imr/(lf+lb) * ( lf*(1-exp(-lb*t1))/lb - lb*(1-exp(-lf*(T-t1)))/lf);

pi2me = k*sig*imr/(lf+lb) * ( (exp(-lb*(t2-T))-exp(-lb*t2))/lb );
x2me =  sig*imr/(lf+lb) * ( lf*(exp(-lb*(t2-T))-exp(-lb*t2))/lb);


pi0mev = sig*imr*onet0*(kv./(lfv+lbv)) .* ...
               (  (exp(-(0-t0)*lfv)-exp(-(T-t0)*lfv)).*(onet0*(1./lfv)));          
x0mev = sig*imr*onet0*(1./(lfv+lbv)) .* ...
               (  -(exp(-(0-t0)*lfv)-exp(-(T-t0)*lfv)).*(onet0*(lbv./lfv)));

pi1mev = sig*imr*onet1*(kv./(lfv+lbv)) .* ...
               ( (1-exp(-t1*lbv)).*(onet1*(1./lbv)) + (1-exp(-(T-t1)*lfv)).*(onet1*(1./lfv)));          
x1mev = sig*imr*onet1*(1./(lfv+lbv)) .* ...
               ( (1-exp(-t1*lbv)).*(onet1*(lfv./lbv)) - (1-exp(-(T-t1)*lfv)).*(onet1*(lbv./lfv)));

pi2mev = sig*imr*onet2*(kv./(lfv+lbv)) .* ...
               ( (exp(-(t2-T)*lbv)-exp(-t2*lbv)).*(onet2*(1./lbv)));          
x2mev = sig*imr*onet2*(1./(lfv+lbv)) .* ...
               ( (exp(-(t2-T)*lbv)-exp(-t2*lbv)).*(onet2*(lfv./lbv)));
           

           
% ---------------------------------------------
% x and pi werning equilibriums, before and after T
% -----------------------------------------

piTw = 0;
Cw = (piTw-k*sig*imr/(lf+lb)*(1-exp(-lb*T))/lb)*exp(lb*T); 

pi0w = Cw*exp(-lb*t0) + pi0me; 
pi1w = Cw*exp(-lb*t1) + pi1me; 
pi2w = Cw*exp(-lb*t2) + pi2me; 

x0w = lf/k*Cw*exp(-lb*t0) + x0me; 
x1w = lf/k*Cw*exp(-lb*t1) + x1me; 
x2w = lf/k*Cw*exp(-lb*t2) + x2me; 

% check explicit formula
if 0; % this is just here as a check on algebra
    pi0wx =  k*sig*imr/(lf+lb)*((exp(lb*(-t0))-exp(lb*(T-t0)))/lb + ...
        (exp(-lf*(-t0))-exp(-lf*(T-t0)))/lf);
    pi1wx =  k*sig*imr/(lf+lb)*((1-exp(lb*(T-t1)))/lb + ...
        (1-exp(-lf*(T-t1)))/lf);
    pi2wx = 0*t2;
    
    x0wx =  k*sig*imr/(lf+lb)*(lf*(exp(lb*(-t0))-exp(lb*(T-t0)))/lb - ...
        lb*(exp(-lf*(-t0))-exp(-lf*(T-t0)))/lf);
    x1wx =  k*sig*imr/(lf+lb)*(lf*(1-exp(lb*(T-t1)))/lb - ...
        lb*(1-exp(-lf*(T-t1)))/lf);
    x2wx = 0*t2;
    
end;

Cwv = (piTw-sig*imr*kv./(lfv+lbv).*(1-exp(-lbv*T))./lbv).*exp(lbv*T); 

pi0wv = (onet0*Cwv).*exp(-t0*lbv) + pi0mev; 
pi1wv = (onet1*Cwv).*exp(-t1*lbv) + pi1mev; 
pi2wv = (onet2*Cwv).*exp(-t2*lbv) + pi2mev; 

x0wv = (onet0*((lfv./kv).*Cwv)).*exp(-t0*lbv) + x0mev; 
x1wv = (onet1*((lfv./kv).*Cwv)).*exp(-t1*lbv) + x1mev; 
x2wv = (onet2*((lfv./kv).*Cwv)).*exp(-t2*lbv) + x2mev; 


% ---------------------------------------------
% no-jump equilibrium 
% -----------------------------------------

piT0 = 0;
C0 = (piT0-k*sig*imr/(lf+lb)*(1-exp(-lf*T))/lf); 

pi10 = C0*exp(-lb*t1) + pi1me; 
pi20 = C0*exp(-lb*t2) + pi2me; 

x10 = lf/k*C0*exp(-lb*t1) + x1me; 
x20 = lf/k*C0*exp(-lb*t2) + x2me; 

% check explicit formulas

if 0 ; 
    pi00x =  k*sig*imr/(lf+lb)*(1-exp(-lf*T))/lf ...
        *(exp(lf*t0)-exp(-lb*t0));
    pi10x =  k*sig*imr/(lf+lb)*((1-exp(-lb*t1))/lb + ...
        (1-exp(-lf*(T-t1)))/lf -(1-exp(-lf*T))/lf*exp(-lb*t1));
    pi20x = k*sig*imr/(lf+lb)*((1-exp(-lb*T))/lb*exp(-lb*(t2-T))- ...
        (1-exp(-lf*T))/lf*exp(-lb*t2));
    
    x10x =  k*sig*imr/(lf+lb)*(lf*(1-exp(-lb*t1))/lb  ...
        -lb*(1-exp(-lf*(T-t1)))/lf -lf*(1-exp(-lf*T))/lf*exp(-lb*t1));
    x20x = k*sig*imr/(lf+lb)*(lf*(1-exp(-lb*T))/lb*exp(-lb*(t2-T))- ...
        lf*(1-exp(-lf*T))/lf*exp(-lb*t2));
    
    
    figure;
    plot([t2;t2], [[ pi10x; pi20x]- [pi10; pi20]]);
    
    figure;
    plot([t2;t2], [[ x10x; x20x]- [x10; x20]]);
    
end;

C0v = (piT0-sig*imr*kv./(lfv+lbv).*(1-exp(-lfv*T))./lfv); 
 
pi10v = (onet1*C0v).*exp(-t1*lbv) + pi1mev; 
pi20v = (onet2*C0v).*exp(-t2*lbv) + pi2mev; 

x10v = (onet1*((lfv./kv).*C0v)).*exp(-t1*lbv) + x1mev; 
x20v = (onet2*((lfv./kv).*C0v)).*exp(-t2*lbv) + x2mev; 


% note the no jump equilibrium does not make any sense as an expected one
% to push the left edge, find new c0 that set no jump at time t-x 

TEv = [-1 -2 -3]; % times for announcement of trap. inflaiton is zero on this date
CEv = exp(lb*TEv).*(0*TEv - sig*imr*k/(lf+lb).*(exp(-lf*(0-TEv))-exp(-lf*(5-TEv)))/lf);

pi0Ev = exp(-lb*t0)*CEv + pi0me*ones(size(CEv));
for j = 1:size(TEv,2); 
    pi0Ev(t0<TEv(j),j)=0*pi0Ev(t0<TEv(j),j);
end; 
pi1Ev = exp(-lb*t1)*CEv + pi1me*ones(size(CEv)); 
pi2Ev = exp(-lb*t2)*CEv + pi2me*ones(size(CEv)); 

x0Ev = lf/k*exp(-lb*t0)*CEv + x0me*ones(size(CEv)); 
for j = 1:size(TEv,2); 
    x0Ev(t0<TEv(j),j)=0*x0Ev(t0<TEv(j),j);
end;
x1Ev = lf/k*exp(-lb*t1)*CEv + x1me*ones(size(CEv)); 
x2Ev = lf/k*exp(-lb*t2)*CEv + x2me*ones(size(CEv)); 


% -----------------------------------------
% All equilibria -- Figure 1
% --------------------------------------------

% to index by inflation at time 0
if 1;
    piT0all = [-3 -0.6 -0.2 -0.05 0 0.05 0.10 0.3 ]; % which eq to show
    Call = (piT0all-k*sig*imr/(lf+lb)*(1-exp(-lf*T))/lf);
end

% index by addition to homogenous solution 
% Call = .... 

onesC = ones(1,size(Call,2));
pi0all = exp(-lb*t0)*Call + pi0me*onesC; 
pi1all = exp(-lb*t1)*Call + pi1me*onesC; 
pi2all = exp(-lb*t2)*Call + pi2me*onesC; 


% -------------------------
% make time series plots
% -------------------------

% --------------------------------------------
% Standard (werning) equilibrium
% --------------------------------------------

figure
plot([-2 10],[0 0], '-k')
hold on
%plot(0,0,'rv','markersize',10,'markerfacecolor','r')
% begin and end with markers not line?
plot([0 5], [5 5], '-k','linewidth',2);
plot([5 5],[-150 10],'--k');
plot([0 0],[-150 10],'--k');
plot(5,0,'rs','markersize',10,'markerfacecolor','r')
plot([t1],100*[pi1wv],'--r','linewidth',1);
plot([-2;0;t1;t2],100*[0; 0; x1w;x2w],'--b','linewidth',3); 
plot([-2; 0; t1;t2],100*[0; 0; pi1w;pi2w],'-r','linewidth',3);

%legend('x','\pi','location','southeast');
ylabel('Percent');
xlabel('Time');
%title('Standard equilibrium, varying \kappa')
text(2,5,' r = -2% trap ','fontsize',16,'backgroundcolor','w')

text(0.2,-20,'\pi -','color','r','fontsize', 20);
text(0.5,-0100,'x--','color','b','fontsize', 18);
text(2.5,-100,'More \kappa \rightarrow','color','r','fontsize',16,'backgroundcolor','w');

plot([5 5],[-150 10],'--k'); %vertical line at trap end
axis([-2 8 -150 10])
print -depsc2 standard_vk.eps
print -dpng standard_vk.png

% --------------------------------------------
% C=0 equilbrium 
% --------------------------------------------

figure
plot([0 10],[0 0], '-k')
hold on
plot([5 5],[-150 10],'--k');
plot([0 0],[-150 10],'--k');
plot(t0,100*x0me,'--b','linewidth',1)
plot(t0,100*pi0me,'-r','linewidth',1)
plot([t0;t1;t2],100*[pi0mev;pi1mev;pi2mev],'--r','linewidth',1);
plot([[-5; 0]; t1;t2],100*[0; 0 ;x1me;x2me],'--b','linewidth',3); 
plot([-5; 0;t1;t2],100*[0; 0;pi1me;pi2me],'-r','linewidth',3);

%plot([5 5],[-1 6],'--k');
axis([-3 8 -1 2.2])
text(0.5,1.1,'\pi','color','r','fontsize', 20);
text(1,-0.5,'x','color','b','fontsize', 18);
text(0.5,2.1,'More \kappa','color','r','fontsize',16);
text(-0.8,0.5,'Expected','color','r','fontsize',16,'horizontalalignment','r')
%legend('x','\pi');
xlabel('Time');
ylabel('Percent');
%title('Backward-stable equilibrium, varying \kappa')
print -depsc2 local_frictionless_vk.eps
print -dpng local_frictionless_vk.png

% --------------------------------------------
% no - jump equilibrium
% --------------------------------------------

figure
plot([-5 10],[0 0], '-k')
hold on
plot([5 5],[-150 10],'--k');
plot([0 0],[-150 10],'--k');
plot([t0; t1; t2], 100*[pi0Ev; pi1Ev; pi2Ev], '-r') 
plot([t0; t1; t2], 100*[x0Ev; x1Ev; x2Ev], '--b') 
plot([[-5;0];t1;t2],100*[0;0;x10;x20],'--b','linewidth',3); 
plot([-5;0;t1;t2],100*[0;0;pi10;pi20],'-r','linewidth',3);
plot(0,0,'rs','markersize',10,'markerfacecolor','r')
plot(-1,0,'rs','markersize',7,'markerfacecolor','r')
plot(-2,0,'rs','markersize',7,'markerfacecolor','r')
plot(-3,0,'rs','markersize',7,'markerfacecolor','r')

text(1,0.7,'\pi','color','r','fontsize',20)
text(1,-1.4,'x','color','b','fontsize',18)
text(-1,0.5,'Expected','color','r','fontsize',16,'horizontalalignment','r')

%plot([t1;t2],100*[pi10v;pi20v],'--r','linewidth',1);

%legend('x','\pi');
ylabel('Percent');
xlabel('Time');
%title('No-jump equilibrium, varying \kappa')
%plot([5 5],[-100 10],'--k');
axis([-3 8 -2 2])
print -depsc2 nojump_vk.eps
print -dpng nojump_vk.png


% -----------------------------------------
% plot of all equilibria 
% --------------------------------------------

figure;
plot([-10 10],[0 0], '-k')
hold on
plot([5 5],[-30 10],'--k');
plot([0 0],[-30 10],'--k');;
plot([t0;t1;t2],100*[pi0all;pi1all;pi2all],'-b');
plot([t1;t2],100*[pi10;pi20],'-r','linewidth',3)
plot([t0;t1;t2],100*[pi0me;pi1me;pi2me],'-r','linewidth',3); % echo the big ones
plot([t0;t1;t2],100*[pi0w;pi1w;pi2w],'-r','linewidth',3);
plot(5,0,'rs','markersize',10,'markerfacecolor','r')
plot(0,0,'rs','markersize',10,'markerfacecolor','r')

axis([-3 8 -15 5])
ylabel('Percent');
xlabel('Time');
%title('Inflation across equilibria')
print -depsc2 inflation_all.eps
print -dpng inflation_all.png


% --------------
% Multipliers
% ------------




%base case multiplier. Not calculated for t<0. 


dCdg0 = -lf/(lf+lb)*(1-exp(-lf*T))*exp(-lb*[t1;t2]); % dC/dg term, pi0=0 solution 
dCdgT = -lf/(lf+lb)*(1-exp(lb*T))*exp(-lb*[t1;t2]); % dC/dg term, piT=0 solution
dx2dg = [ -1 + 1/(lf+lb)*(lf*exp(-lb*t1)+lb*exp(-lf*(T-t1)));...
            -lf/(lf+lb)*(exp(-lb*(t2-T))-exp(-lb*t2))];        % dx2 tg term     

% versions for vector of k %
dCdg0v = -([onet1;onet2]*(lfv./(lfv+lbv).*(1-exp(-lfv*T)))).*exp(-[t1;t2]*lbv); % dC/dg term, pi0=0 solution 
dCdgTv = -([onet1;onet2]*(lfv./(lfv+lbv).*(1-exp(lbv*T)))).*exp(-[t1;t2]*lbv); % dC/dg term, piT=0 solution
dx2dgv =  [ -1 + onet1*(1./(lfv+lbv)) .* ( ((onet1*lfv).*exp(-t1*lbv) + (onet1*lbv).*exp(-(T-t1)*lfv)) );...
               -onet2*(lfv./(lfv+lbv)) .* (exp(-(t2-T)*lbv)-exp(-t2*lbv))];        % dx2 tg term      
        
multf =  dCdg0 + dx2dg; %Multiplier, pi0=0 solution
multw =  dCdgT + dx2dg; %Multipleir, piT=0 solution
multC0= dx2dg;          %Multiplier, C=0 solution

multfv =  dCdg0v + dx2dgv;
multwv =  dCdgTv + dx2dgv;
multC0v= dx2dgv;

figure; 
hold on;
plot([t1;t2],[multw],'-r','linewidth',2); 
plot([t1;t2],[multf],'-b','linewidth',2); 
%plot([t1;t2],[multC0],'-g','linewidth',2); visually so close it's confusing to add 

plot([t1;t2],[multwv],'--r'); 
plot([t1;t2],[multfv],'--b'); 

plot([0 10],[0 0], '-k');
plot([5 5],[-2 25],'--k');
%plot([0 0],[-2 10],'--k');

axis([0 6 -3 15]);
%title('Multiplier');
text(2,7,'Standard \pi_T=0','horizontalalignment','r','fontsize',18,'color','r')
text(3,12,'More \kappa \rightarrow','color','r','fontsize',18,'backgroundcolor','w')
text(1.5,-2.8,'No-jump \pi_0=0','color','b','fontsize',18,'verticalalignment','bottom')
%legend('Standard equilibrium','Backward-stable',2); 
xlabel('Time')
ylabel('Multiplier')
print -depsc2 multiplier.eps


%

% ------------------------
% model with postponed interest rate rise
% ------------------------

%tauv = [ 0 0.5 log(2)/lb 2]' ; % vector of years to postpone interest rate rise

% calcualates tau such that with pi_T+tau = 0, it also happes that pi_0 = 0
object = @(tau_x) (-r_post)*( exp(-lf*T)*(1-exp(-lf*tau_x))/lf-exp(lb*T)*(exp(lb*tau_x)-1)/lb)...
    -imr*((exp(lb*T)-1)/lb-(1-exp(-lf*T))/lf);
tau0 = fzero(object,log(2)/lb); 

tauv = [ 0 0.6 tau0 0.8]' ; % vector of years to postpone interest rate rise



for tauv_index = 1:size(tauv,1);
    tau = tauv(tauv_index);
    
    % must accomodate tau = 0 and tau not 0.01
    taufloor = floor(tau*10)/10; % tau truncated to 0.01 for time index
    t1a = (0:0.01:T)'; % time scale for simulation t<T with postponed rate rise;
    t2a = [(T:0.01:T+taufloor)';T+tau]; % T+tau>t>T last one in case tau not 0.01 divisible
    t3a = [T+tau; (T+taufloor+0.01:0.01:T+5)']; % after it all
      
    Cw = (-sig*imr*k/(lf+lb)*(exp(lb*T)-1)/lb*exp(-lb*(T+tau)) ...
         +sig*r_post*k/(lf+lb)*(exp(lb*tau)-1)/lb*exp(-lb*tau))*exp(lb*(T+tau));
     % sets pi = 0 at T + tau
    C0 = -sig*imr*k/(lf+lb)*(1-exp(-lf*T))/lf ...
         +sig*r_post*k/(lf+lb)*(1-exp(-lf*tau))/lf*exp(-lf*T);
     % sets pi = 0 at t = 0. 
   
     
    pi1ame =  sig*imr*k/(lf+lb) *  ((1-exp(-lb*t1a))/lb + (1-exp(-lf*(T-t1a)))/lf) ...
            - sig*r_post*k/(lf+lb)* (1-exp(-lf*tau))/lf*exp(lf*(t1a-T));
    x1ame =   sig*imr/(lf+lb) *  (lf*(1-exp(-lb*t1a))/lb -lb* (1-exp(-lf*(T-t1a)))/lf) ...
            - sig*r_post/(lf+lb)* (-lb)*(1-exp(-lf*tau))/lf*exp(lf*(t1a-T)); 
    
    pi1a0 = C0*exp(-lb*t1a) + pi1ame;
    pi1aw = Cw*exp(-lb*t1a) + pi1ame;
    x1a0 = lf*C0*exp(-lb*t1a)/k + x1ame;
    x1aw = lf*Cw*exp(-lb*t1a)/k + x1ame;
    
    pi2ame =  sig*imr*k/(lf+lb) *  (exp(lb*T)-1)/lb*exp(-lb*t2a) ...
        - sig*r_post*k/(lf+lb)* ((1-exp(-lb*(t2a-T)))/lb+(1-exp(-lf*(T+tau-t2a)))/lf);
    x2ame =  sig*imr/(lf+lb) *  lf*(exp(lb*T)-1)/lb*exp(-lb*t2a) ...
        - sig*r_post/(lf+lb)* (lf*(1-exp(-lb*(t2a-T)))/lb-lb*(1-exp(-lf*(T+tau-t2a)))/lf);
    
    pi2a0 = C0*exp(-lb*t2a) + pi2ame;
    pi2aw = Cw*exp(-lb*t2a) + pi2ame;
    x2a0 = lf*C0*exp(-lb*t2a)/k + x2ame;
    x2aw = lf*Cw*exp(-lb*t2a)/k + x2ame;
  
    pi3ame =  sig*imr*k/(lf+lb) *  (exp(lb*T)-1)/lb*exp(-lb*t3a) ...
            - sig*r_post*k/(lf+lb)* (exp(lb*tau)-1)/lb*exp(-lb*(t3a-T));    
    x3ame =  sig*imr/(lf+lb) *  lf*(exp(lb*T)-1)/lb*exp(-lb*t3a) ...
            - sig*r_post*k/(lf+lb)*lf* (exp(lb*tau)-1)/lb*exp(-lb*(t3a-T));  
    pi3a0 = C0*exp(-lb*t3a) + pi3ame;
    pi3aw = Cw*exp(-lb*t3a) + pi3ame;
    x3a0 = lf*C0*exp(-lb*t3a)/k + x3ame;
    x3aw = lf*Cw*exp(-lb*t3a)/k + x3ame;  

        xw(:,tauv_index) = [x1aw; x2aw; x3aw];
        piw(:,tauv_index) = [pi1aw; pi2aw; pi3aw];
        
        x0(:,tauv_index) = [x1a0; x2a0; x3a0];
        pi0(:,tauv_index) = [pi1a0; pi2a0; pi3a0];
        
        ta(:,tauv_index) = [t1a;t2a;t3a];
 end
% make time series plot

figure
for tauv_index = 1:4;
    subplot(2,2,tauv_index)
    plot(ta(:,tauv_index),100*x0(:,tauv_index),'--b','linewidth',2);
    hold on
    plot(ta(:,tauv_index),100*pi0(:,tauv_index),'-r','linewidth',2);
    if tauv_index > 1;  % adds old solution. Pretty but too complex
  %      plot(ta(:,1),100*x0(:,1),'--b','linewidth',1);
         plot(ta(:,1),100*pi0(:,1),'-r','linewidth',1);
    end;
 
    plot([0 10],[0 0], '-k')
    plot([5 5],[-3 6],'--k');
    plot([5+tauv(tauv_index) 5+tauv(tauv_index)],[-3 6],'--k');
    plot(0,0,'rs','markersize',10,'markerfacecolor','r')
    if tauv_index == 3; 
        title(['\tau = ' num2str(tauv(tauv_index),'%2.3f')])
        text(2,2,'\pi','color','r','fontsize',18);
        text(1.5,-1,'x','color','b','fontsize',18);
        %text(3,-2,'\pi, \tau=0','color','r','fontsize',18,'horizontalalignment','l');
        text(5.2,2,'\tau','fontsize',18)
    else;
        title(['\tau = ' num2str(tauv(tauv_index),'%2.1f')])
    end;
     axis([0 7 -3 3])
     if tauv_index == 1;
         %legend('x','\pi',3);
     end
     xlabel('Time');
    ylabel('Percent');
end
%title('Local to frictionless solution')
print -depsc2 local_frictionless_mouth.eps
print -dpng local_frictionless_mouth.png

figure
for tauv_index = 1:4;
    subplot(2,2,tauv_index)
    plot(ta(:,tauv_index),100*xw(:,tauv_index),'--b','linewidth',2);
    hold on
    plot(ta(:,tauv_index),100*piw(:,tauv_index),'-r','linewidth',2);
    if tauv_index > 1;  % adds old solution. Pretty but too complex
  %      plot(ta(:,1),100*xw(:,1),'--b','linewidth',1);
        plot(ta(:,1),100*piw(:,1),'-r','linewidth',1);
    end;
        
    xlabel('Time');
    ylabel('Percent');
    plot([0 10],[0 0], '-k')
    plot([5 5],[-50 50],'--k');
    plot([5+tauv(tauv_index) 5+tauv(tauv_index)],[-50 50],'--k');
    plot(T+tauv(tauv_index),0,'rs','markersize',10,'markerfacecolor','r')
   
    %plot([5+tau 5+tau],[-100 10],'--k');
    if tauv_index == 1; 
        axis([0 7 -30 10])
        title(['\tau = ' num2str(tauv(tauv_index),'%2.1f')])
         %legend('x','\pi',4);
    elseif tauv_index == 2;
        axis([0 7 -30 10])
        title(['\tau = ' num2str(tauv(tauv_index),'%2.1f')])
        % legend('x','\pi',4);
    elseif tauv_index == 3;
        title(['\tau = ' num2str(tauv(tauv_index),'%2.3f')])
        text(2,2,'\pi','color','r','fontsize',18);
        text(1,-1.5,'x','color','b','fontsize',18);
        text(3,-2,'\pi, \tau=0','color','r','fontsize',18,'horizontalalignment','l');
        text(5.2,2,'\tau ','fontsize',18)
        axis([0 7 -3 3])
        % legend('x','\pi');
    elseif tauv_index == 4;
         title(['\tau = ' num2str(tauv(tauv_index),'%2.1f')])
        axis([0 7 -10 30])
         %legend('x','\pi',2);
    end; 
end
print -depsc2 standard_mouth.eps
print -dpng standard_mouth.png

