(* ::Package:: *)

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

% Preliminaries
clear
paramS.NoF=4; % Number of fractors
paramS.maxdate=120; % Horizon for rates
nrepeats=3;
nseeds=20;

% forecasts to include
IncCash=true;
IncCashLR=true;
Inc10=false;


% Read in data
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');
DataForecasts10= csvread('G\Input data\Nom ATSM\ATSM inputs\Forecasts_10__S.csv'); % Not released
DataForecasts= csvread('\Input data\Nom ATSM\ATSM inputs\\Cash forecasts - m -extra.csv'); % Not released
DataForecastsLR= csvread('\Input data\Nom ATSM\ATSM inputs\Cash forecasts - m -extra.csv'); % Not released


% Make data matrices
data.Yields=DataYields(2:size(DataYields,1),3:size(DataYields,2));
data.YieldsT=DataYields(1,3:size(DataYields,2));
data.Rf=DataYields(2:size(DataYields,1),2);
data.Returns=DataReturns(2:size(DataReturns,1),2:size(DataReturns,2));
data.Forecasts10=DataForecasts10(2:size(DataForecasts10,1),2:size(DataForecasts10,2));
data.Forecasts10H=DataForecasts10(1,2:size(DataForecasts10,2));
data.Forecasts=DataForecasts(2:size(DataForecasts,1),2:size(DataForecasts,2));
data.ForecastsH=DataForecasts(1,2:size(DataForecasts,2));

[t, q] = size(data.Returns);

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


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

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

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

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


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


% Settings for MLE
% start=vertcat(Mu,vec(param.BVar),vec(param.SigVar));
start=vec(paramS.BVar);
inds_diag=find(eye(paramS.NoF)==1);
inds_nondiag=find(eye (paramS.NoF)~=1);
ub(inds_diag)=1; % upper bounds, somewaht arbitrary
lb(inds_diag)=1e-7;
ub(inds_nondiag)=inf;
lb(inds_nondiag)=-inf;
TypX=start+0.01;

LLF=LL_func _int(start,data,paramS,M, IncCash, Inc10);


% options = optimoptions('fmincon', 'MaxFunEvals',3600,'MaxIter',400,'MaxSQPIter',1e6,'Diagnostics','on','Display','iter','TolFun',1e-3,'TolCon',1e-1,'TolX',1e-3, 'TypicalX', TypX);%, 'DiffMinChange', 1e-10, 'DiffMaxChange', 0.25);
options = optimoptions('fmincon', 'Diagnostics','off','Display','iter','TypicalX', TypX);%, 'DiffMinChange', 1e-10, 'DiffMaxChange', 0.25);

cons=@myconKK_int;
% potentially need extra constraint for stability that eigen are smller 1

% Do the actual minimization starting OLS estimates

x=start;
for i=1:nrepeats
[x,fval,exitflag] = fmincon(@(x0) LL_func _int(x0,data, paramS, M, IncCash, Inc10),x,[],[],[],[],lb,ub,cons,options);
end
xval=x;
fvalmin=fval

% % Do minimisation at random seeds - only do if have surveys
if IncCash==true || Inc10==true
    while  i<nseeds
        Bseed=randn(paramS.NoF,paramS.NoF)*.1;
        for j=1:paramS.NoF
            Bseed(j,j)=randn^2;
        end
        xseed=vec(Bseed);
        [~, ceqseed]=myconKK_int(xseed);
        if ceqseed==0
            for j=1:nrepeats
                [xseed,fvalseed,exitflagseed] = fmincon(@(x0) LL_func _int(x0,data, paramS, M, IncCash, Inc10),xseed,[],[],[],[],lb,ub,cons,options);
                [~, ceq] = myconKK_int(xseed);
                if exitflagseed<0 || ceq > 0
                    fvalseed=1e10
                end
            end
            if fvalseed<fvalmin
                xval=xseed;
                fvalmin=fvalseed
            end
            i=i+1
        end
    end          
end

% Save MLE estimates
paramS.fval=fval;
paramS.BMLE=xval;
paramS.BMLE=reshape(paramS.BMLE,[paramS.NoF,paramS.NoF]);
paramS.MuMLE=mean (M.X)'-paramS.BMLE*mean (M.X)';
% param.SigmaMLE=xval((param.NoF+1)*param.NoF+1:(2*param.NoF+1)*param.NoF,:);
% param.SigmaMLE=reshape(param.SigmaMLE,[param.NoF,param.NoF]);

% Make errors series
Xfit=zeros(t,paramS.NoF);
for i=1 : t
    Xfit(i,:)=paramS.MuMLE'+M.X(i,:)*paramS.BMLE';
end
M.VMLE=M.X(2:t+1,:)-Xfit;

paramS.SigmaMLE=M.VMLE'*M.VMLE/(t-paramS.NoF*2-1)

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

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

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

Bhat=Beta(:,2+paramS.NoF:1+2*paramS.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+paramS.NoF);
Ahat=Beta(:,1);

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

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



% Calcualte An and param.Bn

paramS.Bn=zeros(paramS.maxdate,paramS.NoF);
paramS.Bn(1,:)=-paramS.delta1';
paramS.An=zeros(paramS.maxdate,1);
paramS.An(1,1)=-paramS.delta0;
for i = 2 : paramS.maxdate
    paramS.Bn(i,:)=paramS.Bn(i-1,:)*(paramS.BMLE-paramS.L1)-paramS.delta1';
    paramS.An(i,:)=paramS.An(i-1,:)+paramS.Bn(i-1,:)*(paramS.MuMLE-paramS.L0)+0.5*(paramS.Bn(i-1,:)*paramS.SigmaMLE*paramS.Bn(i-1,:)'+paramS.SigRet)-paramS.delta0;
end

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

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

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

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

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

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


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

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

  save paramM93_S _ 4fac_cash _both _highsd _ 2022 paramS;
%  save M93_S _ 4fac_cash _both _highsd _ 2022 M;
%  clear
% 



