% Matlab program file to accompany
% Christoffersen and Diebold (1997),
% "Cointegration and Long-Horizon Forecasting"
% Main Simulation for Simple example. 
% Calls exgen.m and mlma1.m
% Saves results in exsim.mat

clear;
global dy2 uhat2 Te;

% User-defined parameters
% -----------------------
Ti=20;      % Cut-off for initial values
Te=100;     % Sample size
lam=1;      % Cointegration parameter
sig2e=1;    % Variance of innovation to x(t)
sig2v=1;    % Variance of innovation to y(t)
MC=4;       % Monte Carlo repetitions  
h=50;       % Maximum forecast horizon

% Initialize Matrices
% -------------------
ehatuni1=zeros(MC,h); ehatuni2=zeros(MC,h);
theta=zeros(MC,2);
lamhat=zeros(MC,1);
ehatdv1=zeros(MC,h); ehatdv2=zeros(MC,h);
ehatuv1=zeros(MC,h); ehatuv2=zeros(MC,h);
eghat=zeros(2,h);
ehateg1t=zeros(MC,h); ehateg2t=zeros(MC,h);
dbeta1=zeros(MC,2); dbeta2=zeros(MC,2);
dbetsum=zeros(2,2);

% Monte Carlo Loop
% ----------------
for n=1:MC; if rem(n,100)==0; disp(n); end;
    % Generate the data using the UVAR representation	    
    [y,dy]=exgen(Ti,Te,h,lam,sig2e,sig2v);
    
    % Estimate the Univariate Parameters	
    dy2=dy(:,2);
    theta(n,:)=fmins('mlma1',[-.4 1]); 
   
    % Form Univariate Predictors & Errors
    unihat1=ones(h,1)*y(Te,1);	   
    ehatuni1(n,:)=(y(Te+1:Te+h,1)-unihat1)';    
    unihat2=ones(h,1) * (y(Te,2) + theta(n,1)*uhat2(Te));
    ehatuni2(n,:)=(y(Te+1:Te+h,2)-unihat2)';    

    % UVAR Estimation & Prediction
    x=y(1:Te-1,:);
    beta1=inv(x'*x)*x'*y(2:Te,1);
    beta2=inv(x'*x)*x'*y(2:Te,2);
    beta=[beta1' ; beta2'];
    for i=1:h;
        uvhat(:,i)=beta^i*y(Te,:)';
    end;
    ehatuv=y(Te+1:Te+h,:)-uvhat';
    ehatuv1(n,:)=ehatuv(:,1)';
    ehatuv2(n,:)=ehatuv(:,2)';  

    % DVAR(1) Estimation & Prediction
    dx=dy(1:Te-1,:); 
    dbeta1(n,:)=(inv(dx'*dx)*dx'*dy(2:Te,1))';
    dbeta2(n,:)=(inv(dx'*dx)*dx'*dy(2:Te,2))';
    dbeta=[dbeta1(n,:); dbeta2(n,:)];  
    dbetsum=zeros(2,2);
    for i=1:h;
        dbeti=dbeta^i;
        dbetsum=dbetsum+dbeti;  
        dvhat(:,i)=y(Te,:)' + dbetsum*dy(Te,:)';
    end;     
    ehatdv=y(Te+1:Te+h,:)-dvhat';
    ehatdv1(n,:)=ehatdv(:,1)';
    ehatdv2(n,:)=ehatdv(:,2)';  

    % Cointegration Estimation & Prediction
    egtru1=ones(h,1)*y(Te,1);	
    egtru2=lam*ones(h,1)*y(Te,1);  
    ehateg1t(n,:)=(y(Te+1:Te+h,1)-egtru1)'; 
    ehateg2t(n,:)=(y(Te+1:Te+h,2)-egtru2)'; 

end;

% Compute trace(MSE) across MC reps
% ---------------------------------
truni=std(ehatuni1).^2+std(ehatuni2).^2;
trdv=std(ehatdv1).^2+std(ehatdv2).^2;
truv=std(ehatuv1).^2+std(ehatuv2).^2;
ttreg=std(ehateg1t).^2+std(ehateg2t).^2;


% Compute trace(Tri) across MC reps with TRUE lambda
% --------------------------------------------------
trttuni=zeros(1,h);
trttdv=zeros(1,h);
trttuv=zeros(1,h);
trtteg=zeros(1,h);

trttuni(1)=std(lam*ehatuni1(:,1)-ones(MC,1).*ehatuni2(:,1)).^2+std(ehatuni1(:,1)).^2;
trttdv(1)=std(lam*ehatdv1(:,1)-ones(MC,1).*ehatdv2(:,1)).^2+std(ehatdv1(:,1)).^2;
trttuv(1)=std(lam*ehatuv1(:,1)-ones(MC,1).*ehatuv2(:,1)).^2+std(ehatuv1(:,1)).^2;
trtteg(1)=std(lam*ehateg1t(:,1)-ones(MC,1).*ehateg2t(:,1)).^2+std(ehateg1t(:,1)).^2;

for i=2:h;
    trttuni(i)=std(lam*ehatuni1(:,i)-ones(MC,1).*ehatuni2(:,i)).^2+std(ehatuni1(:,i)-ehatuni1(:,i-1)).^2;
    trttdv(i)=std(lam*ehatdv1(:,i)-ones(MC,1).*ehatdv2(:,i)).^2+std(ehatdv1(:,i)-ehatdv1(:,i-1)).^2;
    trttuv(i)=std(lam*ehatuv1(:,i)-ones(MC,1).*ehatuv2(:,i)).^2+std(ehatuv1(:,i)-ehatuv1(:,i-1)).^2;
    trtteg(i)=std(lam*ehateg1t(:,i)-ones(MC,1).*ehateg2t(:,i)).^2+std(ehateg1t(:,i)-ehateg1t(:,i-1)).^2;  
end;

save exsim truni trdv truv ttreg trttuni trttdv trttuv trtteg sig2e sig2v h lam;