| [bw,b0,aw,a0]=acmub(Y,p,reps,tol1,tol2,model,silent)
|
function [bw,b0,aw,a0]=acmub(Y,p,reps,tol1,tol2,model,silent)
% This function implements approximately median unbiased (MU) estimation of
% the AR(p) model following Andrews, D.W.K and H.-Y. Chen (1994),
% "Approximately Median-Unbiased Estimation of Autoregressive Models,"
% Journal of Business & Economic Statistics, 12(2), 187-204.
%
% No Matlab toolboxes are required.
%
% Features:
% - median unbiased estimation of the sum of autoregressive coefficients
% in the AR(p) model
% - median unbiased model selection
%
% Written by Thomas Maag, maag@mtec.ethz.ch
% VERSION 0.1, APRIL 2009
% THIS IS A PRELIMINARY VERSION, COMMENTS AND QUESTIONS ARE VERY WELCOME.
%
% Input:
% Y Tx1 coumn vector of the time series
% p order of the autoregressive model
% reps number of bootstrap replications at level 2
% tol1 level 1 convergence criterion
% = 0: the model is reestimated until level 2 convergence is
% archieved in the first iteration or 10 rounds have been
% completed (wmax=10)
% > 0: the model is re-estimated for the specified number of
% rounds
% tol2 level 2 convergence criterion: precision in terms of the
% absolute deviation of the median unbiased persistence estimate
% relative to the initial least squares estimator (standard
% value: 0.01)
% model = 0: no constant, no trend
% = 1: constant, no trend
% = 2: constant and linear trend
% silent = 0: progress information is displayed
% = 1: no progress information is displayed
%
% Output:
% bw = [b1 ... bp const trend]': median unbiased ADF coefficient vector
% b0 = [b1 ... bp const trend]': initial OLS ADF coefficient vector
% aw = [a1 ... ap const trend]': median unbiased AR coefficient vector
% a0 = [a1 ... ap const trend]': initial OLS AR coefficient vector
%
% Example:
% [bw,b0,aw,a0]=acmu(Y,3,0,2500,0,0.005,1,0) estimates an AR(3) model
% with a constant and returns the median unbiased estimate of the sum of
% autoregressive coefficients in bw(1,1)
%
wmax = 10;
w = 1;
q = 1;
[b0,a0] = olsar(Y,p,model);
aw = a0;
if silent == 0
disp(['acmu0.1 is computing median unbiased parameter estimates']);
if tol1 > 0, disp(['Number of first level iterations: ' num2str(tol1) ' Number of second level replications: ' num2str(reps)]); end
if tol1 == 0, disp(['Maximum number of first level iterations: ' num2str(wmax) ' Number of second level replications: ' num2str(reps)]); end
disp(['Initial sum of autoregressive coefficients (SAR0): ' num2str(b0(1,1))]);
if model == 0, disp('Initial OLS estimates a0`=[a1 ... ap]:'); end
if model == 1, disp('Initial OLS estimates a0`=[a1 ... ap const]:'); end
if model == 3, disp('Initial OLS estimates a0`=[a1 ... ap const trend]:'); end
disp(a0');
end
while q == 1
if silent == 0, disp(['First level iteration: ' num2str(w)]); end
%Compute median unbiased estimate of b(1,1) given the
%remaining parameters
[bmu,amu,v] = muar(Y,p,aw,a0,reps,tol2,model,silent);
%Compute ols estimates of the remaining parameters given the median
%unbiased estimate
[bw, aw] = murest(Y,p,bmu,model);
aw(1,1) = amu;
if silent == 0
disp(['Updated sum of autoregressive coefficients (SAR0): ' num2str(bw(1,1))]);
disp(['Updated OLS estimates: ']);
disp(aw')
end
w=w+1;
if tol1 == 0
if (v == 2)||(w > wmax), q = 0; end
else
if w > tol1, q = 0; end
end
end
if silent == 0
if abs(bw(1)) < 1, disp(['MU estimates indicate that (trend) stationary model is adequate.']); end
if abs(bw(1)) >= 1, disp(['MU estimates indicate that a unit root model (with drift) should be preferred.']); end
end
end
function [bmu,amu,v]=muar(Y,p,aw,a0,reps,tol,model,silent)
%Generate coefficient vector
beta = flipud(aw(1:p));
sar0 = sum(a0(1:p));
%Initialize values
[T] = size(Y,1);
bsim = NaN(reps,1);
delta = 1;
v=1;
maxv = 10;
s=1;
while (abs(delta)>tol)&&(v < maxv)
Un = randn(T+50,reps);
for j=1:reps
%Generate vector of simulated series
ysim = NaN(T+50,1);
%Set initial values to zero
ysim(1:p,1) = 0;
%Simulate AR process
for i = (p+1):(T+50)
if model == 1
ysim(i,1)=aw(p+1,1)+ysim(i-p:i-1,1)'*beta+Un(i,j);
elseif model == 2
ysim(i,1)=aw(p+1,1)+(i-50)*aw(p+2,1)+ysim(i-p:i-1,1)'*beta+Un(i,j);
else
ysim(i,1)=ysim(i-p:i-1,1)'*beta+Un(i,j);
end
end
%OLS estimation using simulated series
[btemp] = olsar(ysim(51:T+50,1),p,model);
bsim(j)=btemp(1);
end
amu = beta(p);
bmu = sum(beta);
delta = sar0-median(bsim);
if s*delta > 0
beta(p)=beta(p)+delta*abs(s);
else
s=-0.8*s;
beta(p)=beta(p)+delta*abs(s);
end
if silent == 0
disp(['Second level iteration: ' num2str(v)]);
disp(['SAR0: ' num2str(sar0) ' current SAR MU: ' num2str(bmu) ' Delta: ' num2str(delta) ' SAR for next simulation: ' num2str(sum(beta))]);
end
v = v+1;
end
end
function [b,a]=olsar(Y,p,model)
%OLS estimation of the AR(p) model
%b = [constant b1 b2 ... bp]': ADF coefficient vector
%a = [constant a1 a2 ... ap]': AR coefficient vector
%r = [r1 ... rT]': residual vector
[T] = size(Y,1);
X = adfmat(Y,p);
Y = Y(p+1:T,1);
if model == 1
S = [ones(T-p,1)];
elseif model == 2
S = [ones(T-p,1) [1:1:T-p]'];
else
S = [];
end
Z = [X S];
b = inv(Z'*Z)*Z'*Y;
a = adftoar(b,p);
end
function [b,a]=murest(Y,p,bmu,model)
%Generate coefficient vector
[T] = size(Y,1);
X = adfmat(Y,p);
dY = Y-circshift(Y,1).*bmu;
dY = dY(p+1:T,1);
if model == 1
S = [ones(T-p,1)];
elseif model == 2
S = [ones(T-p,1) [1:1:T-p]'];
else
S = [];
end
Z = [X(:,2:end) S];
%Compute OLS estimator b
b = inv(Z'*Z)*Z'*dY;
%Compute AR coefficients
b = [bmu; b];
a = adftoar(b,p);
end
function [a]=adftoar(b,p)
%Computes AR(p) coefficients from ADF coefficients
%b = [b1 b2 ... bp const trend]': ADF coefficient vector
%a = [a1 a2 ... ap const trend]': AR coefficient vector
[k] = size(b,1);
C = -triu(ones(p,p),0);
C(1,:) = 1;
a = [C^(-1)*b(1:p,1); b((p+1):k,1)];
end
function [X]=adfmat(Y,p)
%Computes ADF regressor matrix X for AR(p) model
%X = [Y(t-1) dY(t-1) ... dY(t-p+1)]
[T] = size(Y,1);
A = NaN(T,p);
for i = 1:p
A(:,i) = circshift(Y,i);
end
C = diag(ones(p-1,1),1)-diag(ones(p,1));
C(1,1)=1;
X=A((p+1):T,:)*C;
end |
|
