clear;
startdate=[1971 8];                             %Startdate as [year month] not earlier than [1971 8]
%Code for the Christensen-Diebold-Rudebusch Latent Factor Model
%Creates Figure 3 as output
load zcus;                                      %Zero coupon US yield curve data monthly at maturities 3,6,9,...180 months
m=zcus(12*(startdate(1)-1971)+startdate(2)-7:end,:);
m=m(flipud([size(m,1):-3:1]'),:);               %Go to quarterly data
dt=m(:,1:2);                                    %First two columns are year and month
ns=[1 2 4 8 12 16 20 24 28 32 36 40];
y=m(:,ns+2)';                                   %Select yields at 3, 6 month and 1-10 year maturities
[p,bigt]=size(y);
%------Starting values--------------------------------------%
%Note these assume that certain yields are measured without error
%Just for initialization purposes
lambda=0.5;
for i=1:p;
    Zmat(i,1)=1;
    Zmat(i,2)=(1-exp(-lambda*ns(i)/4))/(lambda*ns(i)/4);
    Zmat(i,3)=Zmat(i,2)-exp(-lambda*ns(i)/4);
end;
sl99=[1 7 12];
alpha=inv(Zmat(sl99,:))*y(sl99,:);  
xmat=[ones(bigt-1,1) alpha(:,1:bigt-1)'];   
bhat=alpha(:,2:bigt)*xmat*inv(xmat'*xmat);
res=alpha(:,2:bigt)-(alpha(:,2:bigt)*xmat*inv(xmat'*xmat)*xmat');
sigs=[sqrt(mean(res(:,1).^2));sqrt(mean(res(:,2).^2));sqrt(mean(res(:,3).^2))];
yfit=Zmat*alpha;
sigsme=sqrt(mean(mean((y-yfit).^2)))*sqrt(12/9);
param0=[reshape(bhat,12,1);sigs;sigsme;lambda];
%------Optimization--------------------------------------------%
options=optimset('Algorithm','active-set');
param=fmincon(@(param)likel(param,y,ns),param0,[],[],[],[],[],[],@(param)mycon(param),options);
%------Filter and smoother-------------------------------------%
mu=param(1:3); phi=reshape(param(4:12),3,3);
Q=diag(param(13:15).^2); H=eye(p)*(param(16)^2); lambda=param(17);
p0g0=reshape(inv(eye(9)-kron(phi,phi))*reshape(Q,9,1),3,3);
a0g0=inv(eye(3)-phi)*mu;
for i=1:p;
    Zmat(i,1)=1;
    Zmat(i,2)=(1-exp(-lambda*ns(i)/4))/(lambda*ns(i)/4);
    Zmat(i,3)=Zmat(i,2)-exp(-lambda*ns(i)/4);
end;
for t=1:bigt;
    if t==1; atgtm1(:,t)=mu+(phi*a0g0); ptgtm1(:,:,t)=(phi*p0g0*phi')+Q; end;
    if t>1; atgtm1(:,t)=mu+(phi*atgt(:,t-1)); ptgtm1(:,:,t)=(phi*ptgt(:,:,t-1)*phi')+Q; end; 
    v=y(:,t)-(Zmat*atgtm1(:,t));
    F=(Zmat*ptgtm1(:,:,t)*Zmat')+H;
    atgt(:,t)=atgtm1(:,t)+(ptgtm1(:,:,t)*Zmat'*inv(F)*v);
    ptgt(:,:,t)=ptgtm1(:,:,t)-(ptgtm1(:,:,t)*Zmat'*inv(F)*Zmat*ptgtm1(:,:,t));
end;

%The smoother
atgbigt=atgt; ptgbigt=ptgt;
for t=bigt-1:-1:1;
    v99=ptgt(:,:,t)*phi'*inv(ptgtm1(:,:,t+1));
    atgbigt(:,t)=atgt(:,t)+(v99*(atgbigt(:,t+1)-atgtm1(:,t+1)));
    ptgbigt(:,:,t)=ptgt(:,:,t)+(v99*(ptgbigt(:,:,t+1)-ptgtm1(:,:,t+1))*v99');
end;

%Construct the P-measure expectation of short rates
for t=1:bigt;
    g(:,1)=atgbigt(:,t);
    for j=2:40; g(:,j)=mu+(phi*g(:,j-1)); end;
    rforec=Zmat(1,:)*g;
    avexp(t,1)=mean(max(0,rforec(1:40)));   %Impose zero bound
end;
tp=[dt y(12,:)' avexp y(12,:)'-avexp];

%Create Figure 3
plot(tp(:,1)+(tp(:,2)/12),tp(:,5),'Linewidth',3);

axis([startdate(1) 2010.5 min(tp(:,5))-0.5 max(tp(:,5))+0.5]);
ylabel('Percentage Points','FontSize',16,'FontName','Times');
set(gca,'FontSize',16,'FontName','Times')