clear all;
load data.txt;
zt=data(:,1);
gt=data(:,2);
T=length(gt);
tt=2;       % tt=1 if constant only; tt=2 if constant and linear trend
boot=20000; % number of bootstrap replications
const=ones(T-1,1);
trend=(1:1:T-1)';
determ=[const trend];
y=[zt(2:T) gt(2:T)];
y0=[zt(1) gt(1)];
x=[determ(:,1:tt) zt(1:T-1) gt(1:T-1)];
beta=inv(x'*x)*(x'*y);
err=y-x*beta;

% Bias estimation by bootstrap
betab1=zeros(boot,2);
betab2=zeros(boot,2);
detb=zeros(boot,2*tt);
for i=1:boot;
    ind=floor(1+(T-1)*rand(T-1,1));
    errb=[zeros(1,2); err(ind,:)];
    errb=[errb(:,1)-mean(errb(:,1)) errb(:,2)-mean(errb(:,2))];
    yyb=[y0; zeros(T-1,2)];
    for j=2:T;
        yyb(j,:)=[determ(j-1,1:tt) yyb(j-1,:)]*beta+errb(j,:);
    end
    yb=yyb(2:T,:);
    xb=[determ(:,1:tt) yyb(1:T-1,:)];
    betab=inv(xb'*xb)*(xb'*yb);
    detb(i,:)=[beta(1:tt,1)' beta(1:tt,2)']-[betab(1:tt,1)' betab(1:tt,2)'];
    betab1(i,:)=[beta(tt+1,1) beta(tt+2,1)]-[betab(tt+1,1) betab(tt+2,1)];
    betab2(i,:)=[beta(tt+1,2) beta(tt+2,2)]-[betab(tt+1,2) betab(tt+2,2)];
end
mbetab1=mean(betab1);
mbetab2=mean(betab2);
mdetb=mean(detb);
bvar1=[beta(tt+1,1) beta(tt+2,1)]+mbetab1; % bias correction for eq. 1
bvar2=[beta(tt+1,2) beta(tt+2,2)]+mbetab2; % bias correction for eq. 2
detvar=[beta(1:tt,1)' beta(1:tt,2)']+mdetb;
disp([bvar1; bvar2])
bvar=[[detvar(1:tt)' detvar(tt+1:2*tt)']; [bvar1(1,1) bvar2(1,1)];...
    [bvar1(1,2) bvar2(1,2)]];
errn=y-x*bvar; % VAR innovations using the bias-corrected estimates
std(errn)

