function [botest] = singleadlbo(y, x, p, d, m)
% Author:Jing Li; Email: lij14@miamioh.edu; July 28 2010.
% Output:the single-equation ADL threshold cointegration test--BO test.
%        See equation (11) of "ADL tests for Threshold Cointegration",
%        Journal of Time Series Analysis, 31 (2010), 241-254
%        for definition of the BO test.
% Input: y--(N times 1) data vector, the dependent variable 
% Input: x--(N times k) data vector, the conditioning variables,
%           which are assumed to be weakly exogenous 
% Input: p--lag number of differenced data, taking integer values 0, 1, 2...
% Input: d--indicator for the determinsitc term. 
%           0 if no intercept term or trend is included;  
%           1 if only intercept term is included; 
%           2 if both intercept term and trend are included.
% Input: m--indicator for model selection. 
%           0 if non-momentum model is used;  
%           1 if momentum model is used;   
% Example: The matlab command a=singleadlbo(y, x, 2, 1, 1) saves the BO
%          testing statistics in a. The dependent variable is y; the conditioning 
%          variables are x; 2 lags of differenced y and differenced x are 
%          included as regressors; the deterministic term includes only 
%          intercept term; the lagged differene of error correction term is used as
%          the threshold variable (momentum model)
% ****Please cite Li and Lee (Journal of Time Series Analysis, 2010) 
%       if you use this code****

ob = length(y); cx = size(x); cx = cx(2);
y1 = [NaN; y(1:end-1)]; x1 = [NaN*ones(1,cx); x(1:end-1,:)];
dy = y - y1; dx = x - x1;
dyp = []; dxp = [];
for j = 1:p,
    dyp = [dyp [NaN*ones(j,1); dy(1:end-j)]];
    dxp = [dxp [NaN*ones(j,cx); dx(1:end-j,:)]];
end

% Estimating the error correction term
gx = [x ones(ob,1)];
e = y - gx*inv(gx'*gx)*(gx'*y);
e1 = [NaN; e(1:(end-1))]; de = e-e1; de1 = [NaN; de(1:(end-1))];

% m=0 for non-momentum model; m=1 for momentum model %
z = e1; 
if m == 1, z = de1; end
se = sort(z);

wbo = []; v = 0.15; jj = round(v*ob); mjj = round((1-v)*ob);
while jj<= mjj,
    i = (z<se(jj));
    y1m = y1.*i; y2m = y1.*(1-i); x1m = []; x2m = [];
    for j = 1:cx, x1m = [x1m x1(:,j).*i]; x2m = [x2m x1(:,j).*(1-i)]; end
    kind = [y1m y2m x1m x2m]; nkind = size(kind); nkind = nkind(2);
    
    % adding the lagged differenced data if p>0 %
    if p==0,data = [dy kind dx];
    else data = [dy kind dx dyp dxp];
    end
    
    % d=1 for intercept term only, d = 2 for intercept term and trend %
    if d==1, data = [data ones(ob,1)];
    elseif d==2, data = [data ones(ob,1) (1:ob)'];
    end
    
    % trimming away missing values
    data(1:(p+1),:) = [];
    
    dep = data(:,1); ind = data(:,(2:end)); 
    ndep = length(dep); nind = size(ind); nind = nind(2);
    iind = inv(ind'*ind);
    b = iind*(ind'*dep);
    re = dep - ind*b;
    sigma = re'*re/(ndep-nind);
    R = eye(nind); R = R(1:nkind,:);
    
    % Wald statistics %
    wbo = [wbo; b'*R'*inv(R*sigma*iind*R')*R*b];     
    
    jj = jj + 1;
end
 
% sup Wald test %
botest=max(wbo);