I'm interested in determining both the slope regression coefficient and plotting regression lines for autocorrelated time-series datasets of rainfall. Specifically, I'd like to identify the best approach in R that would allow me to visualize the regression line on the original (undifferenced) time-series plot when I need to difference the data to remove stationarity (i.e, when d>0 in an arima model).

As a start, I'm exploring the use of auto.arima (from the forecast package) and sarima (from the astsa package) which can output regression coefficients in the presence of autocorrelation.

For example:

  1. Using auto.arima. The 'drift' of -5.009 represents the slope (see http://robjhyndman.com/hyndsight/arima-trends/)

    > min.ar <- auto.arima(dec.yr.mmin$min_prcp)
    > summary(min.ar)
    Series: dec.yr.mmin$min_prcp 
    ARIMA(1,1,0) with drift         
              ar1    drift
          -0.5138  -5.0089
    s.e.   0.2465   5.7986
    sigma^2 estimated as 949:  log likelihood=-57.82
    AIC=121.64   AICc=124.31   BIC=123.34
    Training set error measures:
                         ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
    Training set -0.9479987 28.52129 23.83494 -2.484233 16.12547 0.7957998 -0.2617352
  2. Using sarima to fit the model and output the slope

      > fit.min <- sarima(dec.yr.mmin$min_prcp, 1,1,0,                       reg=dec.yr.mmin$decade)
    initial  value 3.542448 
      iter   2 value 3.488927
      iter   3 value 3.386967
      iter   4 value 3.383464
      iter   5 value 3.382408
      iter   6 value 3.382051
      iter   7 value 3.382024
      iter   8 value 3.382020
      iter   9 value 3.381925
      iter   9 value 3.381925
      iter   9 value 3.381925
      final  value 3.381925 
      initial  value 3.400729 
      iter   2 value 3.399523
      iter   3 value 3.399490
      iter   4 value 3.399488
      iter   4 value 3.399488
      iter   4 value 3.399488
      final  value 3.399488 

3.Output coefficients

      > fit.min$fit$coef
             ar1       xreg 
      -0.5137696 -0.5009045 
  1. For comparison, this is the output from an OLS regression which may give an incorrect slope due to autocorrelation.

      > m3 <- lm(dec.yr.mmin$min_prcp ~ dec.yr.mmin$decade)
      > summary(m3)
      lm(formula = dec.yr.mmin$min_prcp ~ dec.yr.mmin$decade)
          Min      1Q  Median      3Q     Max 
      -45.504  -8.048   1.892  13.650  38.357 
                          Estimate Std. Error t value Pr(>|t|)   
      (Intercept)        1014.1570   319.9461   3.170  0.00807 **
      dec.yr.mmin$decade   -0.4222     0.1580  -2.672  0.02032 * 
      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      Residual standard error: 23.83 on 12 degrees of freedom
      Multiple R-squared:  0.3731,  Adjusted R-squared:  0.3209 
      F-statistic: 7.142 on 1 and 12 DF,  p-value: 0.02032


The output from sarima identifies the slope coefficient and also intercept when d=0. When differencing is required (e.g., ARIMA (1,1,0) as output above), sarima only outputs the slope.

My question: when d = 1 or more, what approaches in R would allow me to add/visualize the regression line onto the original undifferenced time-series plot. Is it possible to derive the fitted values of the regression line, or derive intercept values from sarima/auto.arima or other package?

Many thanks in advance for your suggestions.

  • $\begingroup$ By "regression line", do you mean the in-sample fit? $\endgroup$ – Stephan Kolassa Jul 9 '15 at 14:58
  • $\begingroup$ What's the question exactly? You can use fitted.values(m), where m is the model to extract the in-sample fit to compare to the original data. I'm not sure what you mean by slope and intercept, as a ARIMA(1,1,0) model is not going to be a line. $\endgroup$ – A. Webb Jul 9 '15 at 15:12
  • $\begingroup$ The data are autocorrelated and I'd like to add a regression line onto the time-series plot (I'm not interested in forecasting the data at this stage, just fitting a linear regression). The sarima and auto.arima functions can output the slope coefficient but I'm not sure how best to visualize the line in the time-series plot. There is no intercept from the latter functions when I have to difference the data (d=1) to remove non-stationarity. $\endgroup$ – Richard Jul 9 '15 at 15:22
  • $\begingroup$ Without looking the data I don't know if a deterministic linear trend (which is what I think you are seeking) is the best alternative. For some data a time-varying trend may be more realistic. A non-deterministic trend could be fitted by means of the local-level model (possibly extended with a seasonal component and external regressor variables). You may have a look at this link for further details and examples about what I mean. $\endgroup$ – javlacalle Jul 9 '15 at 21:04
  • $\begingroup$ Thanks for link - will investigate this in a moment. Here are the data I'm currently plotting (copied below). They are mean annual rainfall by decade. I have many more datasets to analyse, but the main issue I'd like to address with the current dataset is how to difference the data and then derive the slope and intercept to allow me to plot on the original undifferenced data. As noted above I'm using arima to derive the slope, but how to find the intercept? $\endgroup$ – Richard Jul 10 '15 at 11:55

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.