(* ::Package:: *)

% This code estiamtes a nominal ATSM using Adrian Crump and Moench (2013) method for
% This version estiamtes the base model and a bootstrapped version
% For Hambur and Haque (2023). Based on code from ACM (2013)


% Preliminaries
param.NoF=4; % number of factors
param.maxdate=120; % futhest yield
nboot=100000; % bootstraps

boot=true; % toggle on or off bootstraping

% Read in data yeilds in levels and monthly returns
DataYields= csvread('\Input data\Nom ATSM\ATSM inputs\Yields_m - extra.csv');
DataReturns= csvread('\Input data\Nom ATSM\ATSM inputs\Returns_m_red - extra.csv');

% Make data matrices
Yields=DataYields(2:size(DataYields,1),3:size(DataYields,2));
Rf=DataYields(2:size(DataYields,1),2);
Returns=DataReturns(2:size(DataReturns,1),2:size(DataReturns,2));

[t, q] = size(Returns);

% Run PCA to extract principal components
[PC, Mean,~,X]= pca1(Yields);

% Take no of factors used in model. First factor is in last row
X=X(:,size(Yields,2)-param.NoF+1:size(Yields,2));
param.PC=PC(:,size(Yields,2)-param.NoF+1:size(Yields,2));
param.Mean=Mean;





% Enforce average positive loading 
temp.neg                        = mean(param.PC) < 0;
X(:, temp.neg)        = -X(:,temp.neg);
param.PC(:, temp.neg) = -param.PC(:,temp.neg);
M.X=X;
clear temp

% Run a VAR of the factors with 1 lag to get variances and residuals
[param.BVar, param.SigVar, M.VVar, param.Mu, ~]=olsvarc(M.X,1);

param.Mu=mean (M.X)'-param.BVar*mean (M.X)'; % almost a small sample adjustment to mean to ensure corect unconditioanl averege short rate


% Do bootstrapping
if boot==true
   [param.BVar,  param.Mu, M.VVar] = Bootstrap(M.X,nboot);
   M.VVar=M.VVar';
end

% LLK of P dynamics
for jj = 1:t
    llk_p(jj,:)=M.VVar(:,jj)'/param.SigVar*M.VVar(:,jj);
end
LLK_P= -t*param.NoF/2*log(2*pi)-1/2*t*log(det(param.SigVar))-1/2*sum(llk_p);



% Regress returns on constant, residuals from VAR and previous factor
Z=[ones(t,1) M.X(1:t,:) M.VVar']; % Regressor series
Z=Z';

Beta=Returns'*Z'/(Z*Z'); % Coefficients

Vret=Returns'-Beta*Z;
SigRet=trace (Vret*Vret')/(q*t); % variance

Bhat=Beta(:,2+param.NoF:1+2*param.NoF); % Beta coeffcients form paper

Bstar=[];
for i=1:q
        Bstar=[Bstar, vec(Bhat (i,:)'*Bhat(i,:))];
end

Bstar=Bstar';

% Other coeffcient groups from paper
Chat=Beta(:,2:1+param.NoF);
Ahat=Beta(:,1);

% Estimate prices of risk
param.L1=(Bhat'*Bhat)\Bhat'*Chat;

param.L0=(Bhat'*Bhat)\Bhat'*(Ahat+0.5*(Bstar*vec(param.SigVar)+SigRet*ones(q,1)));

% Estiamte short rate coefficinets
SrX=[ones(t+1,1) X];
delta=(SrX'*SrX)\SrX'*(Rf./12);
param.delta0=delta(1,:);
param.delta1=delta(2:1+param.NoF,:);

% Calcualte An and param.Bn

param.Bn=zeros(param.maxdate,param.NoF);
param.Bn(1,:)=-param.delta1';
param.An=zeros(param.maxdate,1);
param.An(1,1)=-param.delta0;
for i = 2 : param.maxdate
    param.Bn(i,:)=param.Bn(i-1,:)*(param.BVar-param.L1)-param.delta1';
    param.An(i,:)=param.An(i-1,:)+param.Bn(i-1,:)*(param.Mu-param.L0)+0.5*(param.Bn(i-1,:)*param.SigVar*param.Bn(i-1,:)'+SigRet)-param.delta0;
end


y=zeros(t+1,param.maxdate);
for i=1:t+1
    y(i,:)=param.An'+X(i,:)*param.Bn';
end
M.y=y;

% Calcualte param.Anrf and param.Bn rf
param.Bnrf=zeros(param.maxdate,param.NoF);
param.Bnrf(1,:)=-param.delta1';
param.Anrf=zeros(param.maxdate,1);
param.Anrf(1,1)=-param.delta0;
for i = 2 : param.maxdate
    param.Bnrf(i,:)=param.Bnrf(i-1,:)*(param.BVar)-param.delta1';
    param.Anrf(i,:)=param.Anrf(i-1,:)+param.Bnrf(i-1,:)*(param.Mu)+0.5*(param.Bnrf(i-1,:)*param.SigVar*param.Bnrf(i-1,:)'+SigRet)-param.delta0;
end

yrf=zeros(t+1,param.maxdate);
for i=1:t+1
    yrf(i,:)=param.Anrf'+X(i,:)*param.Bnrf';
end
M.yrf=yrf;


% Calcualte param.Anrf and param.Bn avg exp
param.Bnrfexp=zeros(param.maxdate,param.NoF);
param.Bnrfexp(1,:)=-param.delta1';
param.Anrfexp=zeros(param.maxdate,1);
param.Anrfexp(1,1)=-param.delta0;
for i = 2 : param.maxdate
    param.Bnrfexp(i,:)=param.Bnrfexp(i-1,:)*(param.BVar)-param.delta1';
    param.Anrfexp(i,:)=param.Anrfexp(i-1,:)+param.Bnrfexp(i-1,:)*(param.Mu)-param.delta0;
end

yrfexp=zeros(t+1,param.maxdate);
for i=1:t+1
    yrfexp(i,:)=param.Anrfexp'+X(i,:)*param.Bnrfexp';
end
M.yrfexp=yrfexp;

% Calckaute forward exp short rates
all_per=[1:param.maxdate];
yrf_fwd=fit_short(M.X,param.delta0,param.delta1,param.Mu,param.BVar-eye(param.NoF),all_per,param.NoF,t+1);
M.yrffwd=yrf_fwd;


% LLK of Q dynamics
for jj = 1:t
    llk_q(jj,:)=Vret (:,jj)'/SigRet*Vret(:,jj);
end

LLK_Q= -t*q/2*log(2*pi)-t*q/2*log(SigRet)-1/2*sum(llk_q);
LLK=LLK_Q+LLK_P;
% % Calculate cannonical correlations rank statisitc
 Zcorrel=[ones(t,1) M.X(1:t,:)];
 Zcorrel=Zcorrel';
 Vcoef= (M.VVar*Zcorrel')/(Zcorrel*Zcorrel'); % Regress V on X to get conditional
 VcondX=M.VVar-Vcoef*Zcorrel; % Conditional portions
 Rcoef= (Returns'*Zcorrel')/(Zcorrel*Zcorrel'); % Regress Returns on X to get conditional
 RcondX=Returns'-Rcoef*Zcorrel; % Conditional portions
% 
 QcondX=[VcondX' RcondX'];
 Srank=cov(QcondX);
 Srank11=Srank(1:param.NoF, 1:param.NoF); % Var matrix of VcondX
 Srank21=Srank(param.NoF+1:size(Srank,1),1:param.NoF); % Cov matrix
 Srank22=Srank(param.NoF+1:size(Srank,1),param.NoF+1:size(Srank,1)); % Var matix of RcondX
% 
 rho=eig(inv(Srank11)*Srank21'*inv(Srank22)*Srank21);
 param.rank=-t*log(1-rho (param.NoF,1)^2);
% 
 inv_sig=inv(param.SigVar);
 Var_bet=SigRet*kron(eye(q),inv_sig(param.NoF,param.NoF));
 param.Wald=Bhat (:,param.NoF)'/Var_bet*Bhat(:,param.NoF);

if boot==true
 save paramM93_ 4fac_boot param;
 save M93_ 4fac_boot M;
end
% clear



if boot==false
 save paramM93_ 4fac param;
 save M93_ 4fac M;
end
