2
$\begingroup$

I have physiological time-series data: ~60k observations per channel, ~100 Hz sampling. I will model individual channels with ~20 regressors.

Under OLS, given temporal autocorrelation in the data, standard errors of the beta estimates will be too small (non-iid errors). One method for addressing the iid violation is to model the error term as an ARIMA process.

From readings and prior experience, my understanding is that when fitting a model with ARIMA errors, the betas are typically similar to those obtained from OLS. Given OLS beta estimates are unbiased, this seems like a reasonable expectation. However, for my current data, I find the ARMA-error betas can differ from OLS estimates by large percentages. This appears to be related to how the data are being modeled. (More on that later.)

(a) Without reference to the example data below, under what conditions should I expect the betas from OLS and regression with ARIMA error to meaningfully differ?

I. Below are plots of R example fits for OLS, GLS-REML, "ARMA(0,0) errors", ARMA(1,1) errors, and ARIMA(2,1,5) errors (model selection by auto.arima). For comparison, fit models are superimposed on the response data. The x axis is truncated for clarity.

What is evident is that OLS and ARMA(0,0) yield comparable fits but ARMA(1,1) and ARIMA(2,1,5) differ from OLS and the fits are less convincing. (Code is at the end of the post.)

R plot

II. If I take the fit OLS model and add to it simulated iid Gaussian residuals with the mean and variance of the OLS residuals, fitting a regression model with ARMA(1,1) error now estimates betas very close to the OLS betas. Similarly, if I specify an ARMA(1,1) model but at lags 2, as opposed to at lags 1, the fit is comparable to OLS. (Plots here are from Matlab. The last panel is discussed below.)

Matlab plot

III. For fullness, below are residual plots for OLS fitting. The upper-left plot is the histogram of the standardized residuals. The red line is the observed data. The back line is the theoretical curve.

Matlab plot 2

IV. Given the results from part II, I expect the design matrix may be driving the fit difference between the OLS model and the model with ARMA(1,1)@lags1 error.

The regressors are encoding a pair of lagged response functions. The first is a variable-duration boxcar response function with 10 regressors at delays, $t_i = a+(i-1)2\delta$, for $i=1,2,...,10$ and where $\delta$ is the time between samples. The second response function is a fixed duration, $2\delta$, boxcar with similarly delayed regressors. To avoid multicollinearity, PCA is applied to the design matrix before fitting.

If I respecify the model such that the delays are $t_i = a+(i-1)\delta$, for $i=1,2,...,20$, the fit response, i.e., $\hat{y}$, again resembles the OLS fit model (see above, "full lag model", bottom panel, second figure). Thus, the process of generating the delayed regressors with delays staggered every $2\delta$ appears to be introducing autoregressive heteroscedasticity. (D'oh!)

(b) Is this the correct diagnosis? Would that befoul model fitting for regression with ARIMA errors? If so, how exactly?

(c) One solution is clearly to respecify the model with regressor delays staggered every $\delta$, not $2\delta$. However, is there anything wrong with retaining the current design matrix and only specifying ARIMA terms at nonconsecutive lags greater than 1 (e.g., in Matlab, regARIMA('ARLags',[2,4,6],'MALags',[2,4],...)?

Appendix. Code and model summaries

R code (Matlab yields similar results), part I:

DF <- data.frame(y, B)

fitMod <- lm(y~0+., data = DF)
fitMod.glsREML <- gls(y~0+., data = DF,method="REML")
fitMod.arima00 <- arima(y, xreg = B, include.mean=FALSE, order=c(0,0,0))
fitMod.arima11 <- arima(y, xreg = B, include.mean=FALSE, order=c(1,0,1))
fitMod.arimaA <- auto.arima(y, xreg = B, include.mean=FALSE)

df = data.frame(
1:length(y), 
y, 
B %*% matrix(coefficients(fitMod),ncol=1), 
B %*% matrix(coefficients(fitMod.glsREML),ncol=1), 
B %*% matrix(coefficients(fitMod.arima00),ncol=1),
B %*% matrix(coefficients(fitMod.arima11)[3:length(coefficients(fitMod.arima11))],ncol=1),
B %*% matrix(coefficients(fitMod.arimaA)[6:length(coefficients(fitMod.arimaA))],ncol=1)
)
colnames(df) <- c("t","y","lm","REML","ARMA00","ARMA11","ARMAauto")

png(file="~/blah.png",width=600,height=600)    
gg <- ggplot(df)
gg = gg + geom_line(aes(x=t,y = y),color="#CCCCCC")
gg = gg + geom_line(aes(x=t,y = lm))
gg = gg + scale_x_continuous(limits = c(10000, 12000))
gg = gg + scale_y_continuous(limits = c(-3, 3))
gg = gg + ggtitle("lm")
gg1=gg;
gg <- ggplot(df)
gg = gg + geom_line(aes(x=t,y = y),color="#CCCCCC")
gg = gg + geom_line(aes(x=t,y = REML))
gg = gg + scale_x_continuous(limits = c(10000, 12000))
gg = gg + scale_y_continuous(limits = c(-3, 3))
gg = gg + ggtitle("gls reml")
gg2=gg;
gg <- ggplot(df)
gg = gg + geom_line(aes(x=t,y = y),color="#CCCCCC")
gg = gg + geom_line(aes(x=t,y = ARMA00))
gg = gg + scale_x_continuous(limits = c(10000, 12000))
gg = gg + scale_y_continuous(limits = c(-3, 3))
gg = gg + ggtitle("arma(0,0)")
gg3=gg;
gg <- ggplot(df)
gg = gg + geom_line(aes(x=t,y = y),color="#CCCCCC")
gg = gg + geom_line(aes(x=t,y = ARMA11))
gg = gg + scale_x_continuous(limits = c(10000, 12000))
gg = gg + scale_y_continuous(limits = c(-3, 3))
gg = gg + ggtitle("arma(1,1)")
gg4=gg;
gg <- ggplot(df)
gg = gg + geom_line(aes(x=t,y = y),color="#CCCCCC")
gg = gg + geom_line(aes(x=t,y = ARMAauto))
gg = gg + scale_x_continuous(limits = c(10000, 12000))
gg = gg + scale_y_continuous(limits = c(-3, 3))
gg = gg + ggtitle("auto.arima(2,1,2)")
gg5=gg;
grid.arrange(gg1,gg2,gg3,gg4,gg5, ncol = 1, heights = c(1, 1,1,1,1))
dev.off()

Model summaries, part I:

>     summary(fitMod)

Call:
lm(formula = y ~ 0 + ., data = DF)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.69268 -0.29470  0.05804  0.40363  1.95224 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
stimresp0   0.366478   0.001973 185.713  < 2e-16 ***
stimresp1  -0.409488   0.005429 -75.426  < 2e-16 ***
stimresp2   0.273765   0.011076  24.716  < 2e-16 ***
stimresp3  -0.136721   0.015727  -8.693  < 2e-16 ***
stimresp4  -0.016871   0.020030  -0.842 0.399627    
stimresp5   0.119983   0.020018   5.994 2.06e-09 ***
stimresp6   0.073251   0.020086   3.647 0.000266 ***
stimresp7   0.079315   0.020007   3.964 7.37e-05 ***
stimresp8   0.194574   0.020132   9.665  < 2e-16 ***
stimresp9  -0.019857   0.020057  -0.990 0.322153    
stimresp10  0.036716   0.020089   1.828 0.067603 .  
stimresp11  0.064182   0.020109   3.192 0.001415 ** 
stimresp12  0.044426   0.020246   2.194 0.028213 *  
stimresp13  0.016519   0.020363   0.811 0.417244    
stimresp14  0.003618   0.020587   0.176 0.860503    
stimresp15  0.096590   0.024527   3.938 8.22e-05 ***
stimresp16  0.014653   0.028122   0.521 0.602323    
stimresp17  0.041220   0.030876   1.335 0.181884    
stimresp18 -0.023041   0.033014  -0.698 0.485231    
stimresp19 -0.019120   0.034280  -0.558 0.577011    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5303 on 60717 degrees of freedom
Multiple R-squared:  0.4032,    Adjusted R-squared:  0.403 
F-statistic:  2051 on 20 and 60717 DF,  p-value: < 2.2e-16

>     summary(fitMod.glsREML)
Generalized least squares fit by REML
  Model: y ~ 0 + . 
  Data: DF 
       AIC      BIC    logLik
  95469.44 95658.73 -47713.72

Coefficients:
                Value  Std.Error   t-value p-value
stimresp0   0.3664779 0.00197336 185.71291  0.0000
stimresp1  -0.4094879 0.00542899 -75.42613  0.0000
stimresp2   0.2737650 0.01107631  24.71627  0.0000
stimresp3  -0.1367215 0.01572707  -8.69338  0.0000
stimresp4  -0.0168711 0.02002991  -0.84229  0.3996
stimresp5   0.1199827 0.02001760   5.99386  0.0000
stimresp6   0.0732508 0.02008644   3.64678  0.0003
stimresp7   0.0793154 0.02000742   3.96430  0.0001
stimresp8   0.1945743 0.02013166   9.66509  0.0000
stimresp9  -0.0198572 0.02005673  -0.99005  0.3222
stimresp10  0.0367163 0.02008903   1.82768  0.0676
stimresp11  0.0641823 0.02010905   3.19171  0.0014
stimresp12  0.0444263 0.02024558   2.19437  0.0282
stimresp13  0.0165188 0.02036298   0.81122  0.4172
stimresp14  0.0036178 0.02058679   0.17573  0.8605
stimresp15  0.0965903 0.02452702   3.93812  0.0001
stimresp16  0.0146533 0.02812186   0.52107  0.6023
stimresp17  0.0412198 0.03087645   1.33499  0.1819
stimresp18 -0.0230412 0.03301415  -0.69792  0.4852
stimresp19 -0.0191200 0.03428002  -0.55776  0.5770

 Correlation: 
           stmrs0 stmrs1 stmrs2 stmrs3 stmrs4 stmrs5 stmrs6 stmrs7 stmrs8 stmrs9 stmr10 stmr11
stimresp1  -0.002                                                                             
stimresp2  -0.001  0.001                                                                      
stimresp3  -0.001  0.000  0.002                                                               
stimresp4   0.000  0.000  0.000  0.000                                                        
stimresp5   0.000 -0.001 -0.002  0.000  0.000                                                 
stimresp6   0.000  0.000  0.000  0.000  0.000 -0.001                                          
stimresp7   0.000  0.000  0.000  0.000  0.000  0.001  0.000                                   
stimresp8   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000                            
stimresp9   0.000  0.000 -0.001  0.000  0.000  0.002  0.000  0.000  0.000                     
stimresp10  0.000  0.000 -0.001  0.000 -0.001  0.002  0.000  0.000  0.000  0.000              
stimresp11  0.000 -0.001 -0.002  0.000 -0.001  0.003  0.000  0.000  0.000  0.000  0.000       
stimresp12  0.000 -0.001 -0.001  0.000 -0.001  0.004  0.000  0.000  0.000  0.000  0.001  0.002
stimresp13  0.000  0.000 -0.001  0.000 -0.001  0.004  0.000  0.000  0.000  0.000  0.001  0.002
stimresp14  0.000 -0.001  0.002 -0.001  0.002 -0.006 -0.001  0.000  0.000 -0.001 -0.002 -0.004
stimresp15  0.000  0.003  0.001  0.001  0.000 -0.001  0.000  0.000  0.000  0.000  0.000  0.000
stimresp16  0.000 -0.001 -0.001 -0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
stimresp17  0.000 -0.002  0.002  0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
stimresp18 -0.002 -0.001  0.001 -0.001  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
stimresp19 -0.002 -0.001 -0.001  0.002  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
           stmr12 stmr13 stmr14 stmr15 stmr16 stmr17 stmr18
stimresp1                                                  
stimresp2                                                  
stimresp3                                                  
stimresp4                                                  
stimresp5                                                  
stimresp6                                                  
stimresp7                                                  
stimresp8                                                  
stimresp9                                                  
stimresp10                                                 
stimresp11                                                 
stimresp12                                                 
stimresp13  0.002                                          
stimresp14 -0.002 -0.001                                   
stimresp15  0.000  0.000  0.000                            
stimresp16  0.000 -0.001  0.001 -0.001                     
stimresp17  0.000  0.000 -0.001 -0.001 -0.001              
stimresp18  0.000  0.000  0.001  0.000 -0.003 -0.002       
stimresp19  0.000  0.000  0.001  0.003  0.000  0.000  0.002

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-5.0771729 -0.5556663  0.1094453  0.7610662  3.6810348 

Residual standard error: 0.5303495 
Degrees of freedom: 60737 total; 60717 residual
>     summary(fitMod.arima00)

Call:
arima(x = y, order = c(0, 0, 0), xreg = B, include.mean = FALSE)

Coefficients:
          B1       B2      B3       B4       B5    B6      B7      B8      B9      B10     B11
      0.3665  -0.4095  0.2738  -0.1367  -0.0169  0.12  0.0733  0.0793  0.1946  -0.0199  0.0367
s.e.  0.0020   0.0054  0.0111   0.0157   0.0200  0.02  0.0201  0.0200  0.0201   0.0201  0.0201
         B12     B13     B14     B15     B16     B17     B18     B19      B20
      0.0642  0.0444  0.0165  0.0036  0.0966  0.0147  0.0412  -0.023  -0.0191
s.e.  0.0201  0.0202  0.0204  0.0206  0.0245  0.0281  0.0309   0.033   0.0343

sigma^2 estimated as 0.2812:  log likelihood = -47651.51,  aic = 95345.02

Training set error measures:
                     ME      RMSE       MAE      MPE     MAPE     MASE      ACF1
Training set 0.04913047 0.5302622 0.4209808 73.49066 167.1198 1.447434 0.7558899
>     summary(fitMod.arima11)

Call:
arima(x = y, order = c(1, 0, 1), xreg = B, include.mean = FALSE)

Coefficients:
         ar1     ma1      B1       B2      B3       B4       B5      B6       B7      B8      B9
      0.6884  0.7611  0.1844  -0.0772  0.0652  -0.0304  -0.0278  0.0191  -0.0021  0.0047  0.0351
s.e.  0.0036  0.0022  0.0049   0.0062  0.0076   0.0092   0.0117  0.0178   0.0137  0.0103  0.0127
         B10     B11     B12     B13      B14      B15     B16     B17      B18     B19      B20
      0.0005  0.0205  0.0030  0.0240  -0.0025  -0.0090  0.0241  0.0376  -0.0124  0.0071  -0.0063
s.e.  0.0114  0.0121  0.0117  0.0113   0.0107   0.0144  0.0124  0.0138   0.0149  0.0158   0.0163

sigma^2 estimated as 0.06817:  log likelihood = -4621.8,  aic = 9289.6

Training set error measures:
                  ME      RMSE       MAE       MPE     MAPE      MASE      ACF1
Training set 0.02521 0.2610972 0.1990534 -101.9551 374.7726 0.6843938 0.1676053
>     summary(fitMod.arimaA)
Series: y 
Regression with ARIMA(2,1,2) errors 

Coefficients:
         ar1      ar2      ma1      ma2  drift                                                    
      0.8575  -0.5174  -0.1677  -0.4387  0e+00  0.0692  -0.0618  0.0290  -0.0111  -0.0058  -0.0007
s.e.  0.0058   0.0047   0.0076   0.0069  6e-04  0.0051   0.0054  0.0059   0.0068   0.0091   0.0156

      0.0017  0.0011  0.0128  -0.0106  0.0109  -0.0166  0.0171  -0.0075  0.0001  0.0068  0.0235
s.e.  0.0112  0.0085  0.0105   0.0092  0.0097   0.0093  0.0087   0.0085  0.0118  0.0092  0.0107

      -0.0101  0.0098  -0.0101
s.e.   0.0121  0.0135   0.0144

sigma^2 estimated as 0.0589:  log likelihood=-169.75
AIC=391.5   AICc=391.53   BIC=625.88

Training set error measures:
                        ME      RMSE       MAE       MPE     MAPE      MASE       ACF1
Training set -1.636533e-05 0.2426426 0.1769946 -77.68525 313.9655 0.6085504 0.05154149
> 
> 

Matlab code, part II:

OLSstatsX = fitlm(X, y, [eye(size(X,2)) zeros(size(X,2),1)]);

fh = figure('visible',visibility);
set(fh,'Color','w');
set(fh,'Units','normalized');
set(fh,'position',[screensize(1:2) screensize(3)/2 screensize(4)/2]);

[EstMdlX11,EstParamCov,logL1,info] = estimate(regARIMA('AR',0,'MA',0,'intercept',0),y,'X',X,'Display','params','beta0',OLSstatsX.Coefficients{:,1});

subplot(5,1,1);
plot(y(10000:12000),'color',[0.7 0.7 0.7])
hold on
plot(X((10000:12000),1:(nBase))*EstMdlX11.Beta(:),'k')
set(gca,'ylim',[-3 3])
title('real data - ARMA(0,0)')
set(gca,'xlim',[0 2000])

[EstMdlX11,EstParamCov,logL1,info] = estimate(regARIMA('ARLags',[1],'MALags',[1],'intercept',0),y,'X',X,'Display','params','beta0',OLSstatsX.Coefficients{:,1});

subplot(5,1,2);
plot(y(10000:12000),'color',[0.7 0.7 0.7])
hold on
plot(X((10000:12000),1:(nBase))*EstMdlX11.Beta(:),'k')
set(gca,'ylim',[-3 3])
title('real data - ARMA(1,1) at lag 1')
set(gca,'xlim',[0 2000])

e = OLSstatsX.Residuals.Raw;    
z = X*OLSstatsX.Coefficients{:,1}+normrnd(mean(e),std(e),size(e));
[EstMdlX11,EstParamCov,logL1,info] = estimate(regARIMA('ARLags',[1],'MALags',[1],'intercept',0),z,'X',X,'Display','params','beta0',OLSstatsX.Coefficients{:,1});

subplot(5,1,3);
plot(z(10000:12000),'color',[0.7 0.7 0.7])
hold on
plot(X((10000:12000),1:(nBase))*EstMdlX11.Beta(:),'k')
set(gca,'ylim',[-3 3])
title('fake data - ARMA(1,1) at lag 1')
set(gca,'xlim',[0 2000])

[EstMdlX11,EstParamCov,logL1,info] = estimate(regARIMA('ARLags',[2],'MALags',[2],'intercept',0),y,'X',X,'Display','params','beta0',OLSstatsX.Coefficients{:,1});

subplot(5,1,4);
plot(y(10000:12000),'color',[0.7 0.7 0.7])
hold on
plot(X((10000:12000),1:(nBase))*EstMdlX11.Beta(:),'k')
set(gca,'ylim',[-3 3])
title('real data - ARMA(1,1) at lag 2')
set(gca,'xlim',[0 2000])

Xc = [ Xb [Xb(end,:) ; Xb(1:end-1,:)] ];
[~,C,~,~,~,~] = pca(Xc,'centered',false);
C(isnan(C)) = 0;
C = C(include_t,:);
[EstMdlX11,EstParamCov,logL1,info] = estimate(regARIMA('ARLags',[1],'MALags',[1],'intercept',0),y,'X',C,'Display','params');

subplot(5,1,5);
plot(y(10000:12000),'color',[0.7 0.7 0.7])
hold on
plot(C((10000:12000),:)*EstMdlX11.Beta(:),'k')
set(gca,'ylim',[-3 3])
title('full lag model, real data - ARMA(1,1) at lag 1')
set(gca,'xlim',[0 2000])
$\endgroup$
  • $\begingroup$ This is getting slightly off topic so I am adding it as a comment but, bonus question! (d) Alternately one can use HAC estimates of the standard error. If the option presented in (c) were inadvisable, i.e., ARMA(1,1) errors at lags 2 with the current design matrix, would HAC standard errors be similarly inadvisable? $\endgroup$ – jbjo Aug 22 '18 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.