%This code is the likelihood function to be minmised to estiamte the risk free dynamics  parameters for a 
%nominal terms structure model using Adrian Crump and Moench (2013) method but including surveys.
% Jonathan Hambur 2016

function W = LL_func(x0,data, param,M, IncCash, Inc10)

[t, f] = size(M.X);
[~, F10] = size(data.Forecasts10);
[~, F] = size(data.Forecasts);
[~, q] = size(data.Returns);

%Make the matricies that we use later


 %coefficeints on factors
Coeff=reshape(x0(param.NoF+2:end),[param.NoF,param.NoF]);
Mu=mean(M.X)'-param.BVar*mean(M.X)';

r0=x0(1);
r1=x0(2:param.NoF+1);

T=param.maxdate;

Ret=data.Returns;

% Sigma=x0(param.NoF*(param.NoF+1)+1:(param.NoF)*(param.NoF*2+1));
% Sigma=reshape(Sigma,[param.NoF,param.NoF]);

%Set the errors 
%sd's of bond rates forecasts and excess returns - 10 bps and 50 bps
forc_sd = 0.005;

%Make inputs cash
Th=data.ForecastsH;
Th_p=round(data.ForecastsH*12);
For=data.Forecasts./100;

%Make inputs 10-year
Th10=data.Forecasts10H;
Th_p10=round(data.Forecasts10H*12);
For10=data.Forecasts10./100;


% LL terms
LLX=0;
LLY=0;
LLY10=0;
LLRF=0;

%% Log-liklihood function for PCs
%Calculate fitted PC and errors
Xfit=zeros(t-1,f);
for i= 1:t-1
    Xfit(i,:)=Mu'+M.X(i,:)*Coeff';
end
Xerr=M.X(2:end,:)-Xfit;
Xerr=Xerr-ones(t-1,1)*mean(Xerr);

Sigma = Xerr'*Xerr/(t-1-f*2);

%Scaled X erros
Xscale=zeros(t-1,1);
for i = 1:t-1
    Xscale(i)=Xerr(i,:)/Sigma*Xerr(i,:)';
end

% Likelihood - ignore pi term as constant
LLX=  (t-1)/2*log(det(Sigma)) + .5.*sum(Xscale); % 1*T



%% Log-likelihood function for the 10-year bonds surveys
 
%%%%% Estiamte other parameters using least squares

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

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

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

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

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

L0=(Bhat'*Bhat)\Bhat'*(Ahat+0.5*(Bstar*vec(Sigma)+ret_var*ones(q,1)));

if Inc10==true
    
    % % % Make matrices of iterated coeffcients for bond pricing measure with risk
    [H0, H1] = Iterate(Mu,Coeff,Sigma,r0, r1,T,f,ret_var, L0, L1); %Fuction that computes A and B matrices for maturities
    Y=10; %Forecast for Y year rate
    k=Y*12; % In quarters
    yfit10=zeros(t,F10);
    yerr10=zeros(t,F10);
    ysum10=zeros(1,F10); %matrix to sum errors

    all_per10 = [1:max(Th_p10)+1];

    y_all10= fit_short(M.X,H0(k),H1(k,:)',Mu,Coeff-eye(param.NoF),all_per10,param.NoF,t);

    for i=1:t
        for j=1:F10;
            if For10(i,j)==0
            else
                yfit10(i,j)=-y_all10(i,Th_p10(j)+1)/10;
                yerr10(i,j)=For10(i,j)-yfit10(i,j);
                ysum10(1,j)=ysum10(1,j)+(yerr10(i,j)).^2/(forc_sd^2*Th10(j));
            end
        end
    end

    %Log-liklihood function (no need to put in variance of pi term as constant
    LLY10=0.5*sum(ysum10');
end

%% Log-likelihood function for the cash rate (taken as 1-motnh risk-free)
% Make matrices of iterated coeffcients for bond pricing measure without risk
% [H0f, H1f] = Iteratef(Mu,Coeff,Sigma,r0, r1,T,f,ret_var); %Fuction that computes A and B matrices for maturities

if IncCash==true;
    Y=1/12; %Forecast for Y year rate
    k=Y*12; % In months
    yfit=zeros(t,F);
    yerr=zeros(t,F);
    ysum=zeros(1,F); %matrix to sum errors

    all_per = [1:max(Th_p)+1];

    y_all= fit_short(M.X,r0,r1,Mu,Coeff-eye(param.NoF),all_per,param.NoF,t);

    for i=1:t
        for j=1:F;
            if For(i,j)==0
            else
                yfit(i,j)=y_all(i,Th_p(j)+1)*12;
                yerr(i,j)=For(i,j)-yfit(i,j);
                ysum(1,j)=ysum(1,j)+(yerr(i,j)).^2/(forc_sd^2*Th(j));
            end
        end
    end


    %Log-liklihood function (no need to put in variance of pi term as constant
    LLY=0.5*sum(ysum');
end


%% Log-liklihood function for RF rates
%Calculate fitted rf rates and errors
rffit=zeros(t,1);
for i= 1:t
    rffit(i,:)=r0+M.X(i,:)*r1;
end
rferr=data.Rf/12-rffit;
rferr_dem=rferr-mean(rferr);

Sigmarf = rferr_dem'*rferr_dem/(t-1);



% Likelihood - ignore pi term as constant
LLRF=  (t)/2*log(Sigmarf) + .5/Sigmarf.*sum(rferr.^2); % 1*T
%% Overall LL 
W=LLX+LLY+LLY10+LLRF;
            
            
            
            
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [An, Bn]=Iterate(Mu,Coeff,Sigma,r0, r1,T,f,ret_var, L0, L1)
Bn=zeros(T,f);
Bn(1,:)=-r1';
An=zeros(T,1);
An(1,:)=-r0;
for i= 2 : T
    Bn(i,:)=Bn(i-1,:)*(Coeff-L1)-r1';
    An(i,:)=An(i-1,:)+Bn(i-1,:)*(Mu-L0)+0.5*(Bn(i-1,:)*Sigma*Bn(i-1,:)'+ret_var)-r0;
end

% function [An, Bn]=Iteratef(Mu,Coeff,Sigma,r0, r1,T,f,ret_var)
% Bn=zeros(T,f);
% Bn(1,:)=-r1';
% An=zeros(T,1);
% An(1,:)=-r0;
% for i= 2 : T
%     Bn(i,:)=Bn(i-1,:)*(Coeff)-r1';
%     An(i,:)=An(i-1,:)+Bn(i-1,:)*(Mu)+0.5*(Bn(i-1,:)*Sigma*Bn(i-1,:)'+ret_var)-r0;
% end
