matlab中用AR项估计线性回归的快速方法

Posted 2023-03-25

【中文标题】matlab中用AR项估计线性回归的快速方法

【英文标题】：Fast approach in matlab to estimate linear regression with AR terms

【发布时间】：2020-08-08 05:49:21

【问题描述】：

我正在尝试使用 AR 误差项估计（负载）线性回归的回归和 AR 参数。 （您也可以将其视为具有外生变量的 MA 过程）：

, 其中 , 滞后长度为 p

我正在关注official matlab recommendations 并使用regArima 设置一些回归并提取回归和AR 参数（参见下面的可重现示例）。

问题： regArima 很慢！对于 5 次回归，matlab 需要14.24sec。我打算运行大量不同的回归模型。有没有更快的方法？

y = rand(100,1);
r2 = rand(100,1);
r3 = rand(100,1);
r4 = rand(100,1);
r5 = rand(100,1);
exo = [r2 r3 r4 r5];

tic
 for p = 0:4
     Mdl = regARIMA(3,0,0);
     [EstMdl, ~, LogL] = estimate(Mdl,y,'X',exo,'Display','off');

 end
toc

【参考方案1】：

与使用最大似然法的 regArima 函数不同，Cochrane-Orcutt 程序依赖于 OLS 回归的迭代。当这种方法有效时，还有一些特殊性（请参阅发布的链接）。但是对于这个问题的目的，方法是有效的，而且很快！

我修改了 James Le Sa​​ge 的代码，该代码仅涵盖 1 阶的 AR 滞后，以涵盖 p 阶的滞后。

function result = olsc(y,x,arterms) 
% PURPOSE: computes Cochrane-Orcutt ols Regression for AR1 errors 
%--------------------------------------------------- 
% USAGE: results = olsc(y,x) 
% where: y = dependent variable vector (nobs x 1) 
%        x = independent variables matrix (nobs x nvar) 
%--------------------------------------------------- 
% RETURNS: a structure 
%        results.meth  = 'olsc' 
%        results.beta  = bhat estimates 
%        results.rho   = rho estimate 
%        results.tstat = t-stats 
%        results.trho  = t-statistic for rho estimate 
%        results.yhat  = yhat 
%        results.resid = residuals 
%        results.sige  = e'*e/(n-k) 
%        results.rsqr  = rsquared 
%        results.rbar  = rbar-squared 
%        results.iter  = niter x 3 matrix of [rho converg iteration#] 
%        results.nobs  = nobs 
%        results.nvar  = nvars 
%        results.y     = y data vector 
% -------------------------------------------------- 
% SEE ALSO: prt_reg(results), plt_reg(results) 
%--------------------------------------------------- 

% written by: 
% James P. LeSage, Dept of Economics 
% University of Toledo 
% 2801 W. Bancroft St, 
% Toledo, OH 43606 
% [e-mail] 

% do error checking on inputs 
if (nargin ~= 3); error('Wrong # of arguments to olsc'); end; 

[nobs nvar] = size(x); 
[nobs2 junk] = size(y); 

if (nobs ~= nobs2); error('x and y must have same # obs in olsc'); end; 

% ----- setup parameters 
ITERMAX = 100; 
converg = 1.0; 
rho = zeros(arterms,1); 
iter = 1; 
% xtmp = lag(x,1); 
% ytmp = lag(y,1); 

% truncate 1st observation to feed the lag 
% xlag = x(1:nobs-1,:); 
% ylag = y(1:nobs-1,1); 
yt = y(1+arterms:nobs,1); 
xt = x(1+arterms:nobs,:); 

xlag = zeros(nobs-arterms,arterms);
for tt = 1 : arterms
xlag(:,nvar*(tt-1)+1:nvar*(tt-1)+nvar) = x(arterms-tt+1:nobs-tt,:);
end

ylag = zeros(nobs-arterms,arterms);
for tt = 1 : arterms
ylag(:,tt) = y(arterms-tt+1:nobs-tt,:);
end

% setup storage for iteration results 
iterout = zeros(ITERMAX,3); 

while (converg > 0.0001) & (iter < ITERMAX), 
% step 1, using intial rho = 0, do OLS to get bhat 
 ystar = yt - ylag*rho; 
 xstar = zeros(nobs-arterms,nvar);
 for ii = 1 : nvar
    tmp = zeros(1,arterms);
    for tt = 1:arterms
        tmp(1,tt)=ii+nvar*(tt-1);
    end       
    xstar(:,ii) = xt(:,ii) - xlag(:,tmp)*rho; 
 end

beta = (xstar'*xstar)\xstar' * ystar; 

e = y - x*beta; 

% truncate 1st observation to account for the lag 
et = e(1+arterms:nobs,1);
elagt = zeros(nobs-arterms,arterms);
for tt = 1 : arterms
    elagt(:,tt) = e(arterms-tt+1:nobs-tt,:);
end

% step 2, update estimate of rho using residuals 
%         from step 1 

res_rho = (elagt'*elagt)\elagt' * et;  
rho_last = rho; 
rho = res_rho; 
converg = sum(abs(rho - rho_last)); 

% iterout(iter,1) = rho; 
iterout(iter,2) = converg; 
iterout(iter,3) = iter; 

iter = iter + 1; 

end; % end of while loop 

if iter == ITERMAX 
%  error('ols_corc did not converge in 100 iterations'); 
  print('ols_corc did not converge in 100 iterations'); 

end; 

result.iter= iterout(1:iter-1,:); 

% after convergence produce a final set of estimates using rho-value 
 ystar = yt - ylag*rho; 
 xstar = zeros(nobs-arterms,nvar);
 for ii = 1 : nvar
    tmp = zeros(1,arterms);
    for tt = 1:arterms
        tmp(1,tt)=ii+nvar*(tt-1);
    end       
    xstar(:,ii) = xt(:,ii) - xlag(:,tmp)*rho; 
 end

result.beta = (xstar'*xstar)\xstar' * ystar; 
e = y - x*result.beta; 

et = e(1+arterms:nobs,1);
elagt = zeros(nobs-arterms,arterms);
for tt = 1 : arterms
    elagt(:,tt) = e(arterms-tt+1:nobs-tt,:);
end

u = et - elagt*rho;
result.vare = std(u)^2;

result.meth = 'olsc'; 
result.rho = rho; 
result.iter = iterout(1:iter-1,:); 
% % compute t-statistic for rho 
% varrho = (1-rho*rho)/(nobs-2); 
% result.trho = rho/sqrt(varrho); 

（我没有在最后 2 行中调整长度为 p 的 rho 向量的 t 检验，但这应该是直截了当的......）