function [e,k1,vrat,arcs] = whit1(y,nhi,k2)
% whit1: fit AR model to a time series, returning model information and residuals
% [e,k1,vrat,arcs] = whit1(y,nhi,k2);
% Last revised 7-12-99
%
% You specify the highest order AR model to consider. Models up to that order are
% fit, and the Akaike Information Criterion (AIC) is used to pick the best model.
%
%*** IN **********************
%
% y (my x 1)r time series, vector; NaNs not allowed
% nhi (1 x 1)i highest AR order to consider (if k2==1), or the only
% AR order model to try fitting (if k2=2)
% k2 (1 x 1)i option for order selections
% ==1 fit models of order 1 to nhi
% ==2 fit model of order nhi only
%
%*** OUT **************************
%
% e (my x 1)r AR residuals, with mean added back in; or, if null model, the
% original series y
% k1 - the order of the AR model deemed best by the AIC, or 0 if null model.
% vrat - the ratio of variance of AR residuals to variance of original time series y
% arcs the ar coefs and their two-standard errors in a two-row array; row 1 has
% the coefficients; row two has 2*standard error of the coefficients
%
%*** REFERENCES
%
% The loss function and AIC are discussed in Ljung, L. 1995. System identification
% toolbox; for use with MATLAB, The MathWorks, Inc., p. 3-46.
%
%*** UW FUNCTIONS CALLED -- none
%*** TOOLBOXES NEEDED -- system identification
%
%*** NOTES
%
% vrat: the variance ratio is computed over years k1+1 to n1, where n1 is length of y
%
% null model (k1==0):
% The "Loss function" V is the normalized sos of prediction errors, or the
% variance of the AR residuals.
% A "null model" loss function is defined as the mean sos of the
% variable z (original data with mean subtracted) over the years nhi+1 to the
% last year of z. If the null-model loss function is no higher than the V
% for any estimated model, a "null model" is assumed. This means no AR model is
% justified, and k1==1 is returned.
%
mn=mean(y);
z=dtrend(y);
n1=length(y);
if k2==1; % mode 1, fit up to order nhi
NN=(1:nhi)'; % orders to try
else; % mode 2, fit order nhi only
NN=nhi;
end
V=arxstruc(z,z,NN);
[nn,vmod]=selstruc(V,'AIC');
th=arx(z,nn);
[e,r]=resid(z,th);
k1=nn(1);
nsub=n1-nhi;
vmodnull = log(((nsub-1)/nsub) * (std(z(nhi+1:n1))) ^ 2);
vrat= (std(e(k1+1:n1)) ^2) / (std(z(k1+1:n1)) ^2);
e=e+mn;
thsub=diag(th(4:3+k1,1:k1)); % variances of est params
stderr= 2*sqrt(thsub);
arcs=[th(3,1:k1); stderr'];
% Check that null model (ar order =0) is not better than any of others
% considered.
L1=vmodnull(:,ones(length(vmod),1)) <= vmod(1,:);
if all(L1); % Loss function for null model is smaller than for
% any other model, so do not whiten. Pass the orig
% series back as the "whitened" series, and tell
% that model is null
k1=0;
e=y;
vrat=1.0;
arcs=[0 0]';
end