%=========================================================================%
% This program calculates error detection probabilities building from
% Hansen and Sargent (2008)
%=========================================================================%
% A is the vector of coefficients in the DGP
% MU is the vector of drifts 
% SIGMA is the vector of volatility
% ETHA is the bound for the worst case shock Delta (D)
% T is the dimension of the time series  
% nsim is the number of simulations to perform
%=========================================================================%
function [pa, pb, detecp]=EDP1(A,SIGMA,T,MU,ETHA,nsim)


% Initialize matrices for Monte Carlo simulations
x=zeros(nsim+1,T);
x2=zeros(nsim+1,T);
v=zeros(nsim+1,T);
v2=zeros(nsim+1,T);
cross_vw=zeros(nsim+1,T);
cross_veps=zeros(nsim+1,T);
r_A=zeros(1,nsim+1); % likelihood ratios
r_B=zeros(1,nsim+1); % likelihood ratios 

%=========================================================================%
% PART ONE - Simulate nsim drifts given ETHA
%=========================================================================%
hatETHA=ETHA*rand(1,nsim);
% UDMU=MU+SIGMA*sqrt(2*hatETHA); can be used if want optimistic case
LDMU=MU-SIGMA*sqrt(2*hatETHA);
D=LDMU/MU-1;
MU2=min(LDMU);
D=[0 D];
 
%=========================================================================%
% PART TWO - Generate returns assuming the approximating model is right
%=========================================================================%
for ii=1:nsim+1 
  x(ii,1)=MU;
  x2(ii,1)=MU2;
  wn=0.4*randn(1,T);  %adjusted to make it in the range [-1,1]
  epsn=0.4*randn(1,T); %adjusted to make it in the range [-1,1]
  
  % Generate time series of size T 
  for tt=1:T-1
    x(ii,tt+1)=(ones-A)*MU+A*x(ii,tt)+SIGMA*wn(:,tt);
  end;

  % Calculate likelihood ratio r_A
  v(ii,:)=D(:,ii)*x(ii,:);
  tmp1=v(ii,:).^2;
  tmp2=sum(tmp1,1);
  tmp3=mean(tmp2);
  
  pt11=tmp3/2;
  
  for tt=1:T-1
      cross_vw(ii,tt)=v(ii,tt)'*wn(:,tt);
  end;
  
  pt12=mean(cross_vw(ii,:));
  
  r_A(:,ii)=pt11-pt12;

  %=======================================================================%
  % PART THREE - Generate returns assuming the worst case model is right
  %=======================================================================%
  % Generate time series of size T
  % Compute v using worst case drift for given all t / 1st part of r_A
  % Compute cross vw for given all t / 2nd part of r_A
  for tt=1:T-1
    x2(ii,tt+1)=(ones-A)*MU2+A*x2(ii,tt)+SIGMA*epsn(:,tt);
  end
  
  % Calculate likelihood ratio r_B
  v2(ii,:)=D(:,ii)*x2(ii,:);
  tmp4=v2(ii,:).^2;
  tmp5=sum(tmp4,1);
  tmp6=mean(tmp5);
  
  pt21=tmp6/2;

  for tt=1:T-1
      cross_veps(ii,tt)=v2(ii,tt)'*epsn(:,tt);
  end;
  
  pt22=mean(cross_veps(ii,:));

  r_B(:,ii)=pt21+pt22;

end;

%=========================================================================%
% PART FOUR - Calculating probabilities of making errors type I and II
%=========================================================================%
% Calculating probability of making error I - when A is true
na=0;
nb=0;
for jj=1:nsim
  if r_A(jj) <= 0
    na=na+1;
  end
end

pa=na/nsim;

% Calculating probability of making error II - when B is true
for jj=1:nsim
  if r_B(jj) < 0
    nb=nb+1;
  end
end

pb=nb/nsim;

% Error detection probability under complete ignorance
detecp=(1/2)*(pa+pb);