REG_ARIMA

   more...

Fits a univariate, non‑seasonal ARIMA time series model with the inclusion of one or more regression variables.

Required Arguments

Y — Array of length NOBS containing the time series. (Input)

IMODEL — Array of length 3 containing the model order parameters. (Input)

 

IMODEL(I)	Description
1	Order of the autoregressive part, p, where p≥0.
2	Order of the non‑seasonal difference operator, d, where d≥0.
3	Order of the moving average part, q, where q≥0.

If p = 0 and q = 0, only regression is performed.

PARMA —Array of length 1 + p + q containing the estimated autoregressive (AR) and moving average (MA) parameters of the ARIMA(p,d,q) model. PARMA(1) is the estimated AR constant parameter, PARMA(2: (p+1)) contains the AR parameter estimates and PARMA((p+2):), contains the MA parameter estimates.   (Output)

Optional Arguments

NOBS — Number of observations. (Input)
Default: NOBS = size(Y).

X — Array of size NOBS by K containing the regression data, where K = size(X,2) is the number of user supplied regression variables. (Input)
Specific columns in X may be selected using the INDX argument. Otherwise, all columns of X are used.
Default: No regression variables are included.

XLEAD — Array of size MXLEAD by K containing the regression data to be used in obtaining forecasts, where K = size(X,2) is the number of user supplied regression variables. (Input)
Specific columns in XLEAD may be selected using the INDX argument. Otherwise all columns of XLEAD are used.
Note: If MXLEAD > 0 and optional argument X is present, XLEAD is required.

INDX — Index array containing the column numbers in X and XLEAD that are to be used for the regression variables. (Input)
Default: All columns of X and XLEAD are used.

TREND — Logical. If .TRUE., the routine will include a trend variable. (Input)
Note: Setting TREND = .TRUE. has the effect of fitting an intercept term in the regression. If the difference operator IMODEL(2) = d > 0, the effect on the model in the original, undifferenced space is polynomial of order d.
Default: TREND = .TRUE..

MXLEAD — Maximum lead time for forecasts. (Input)
Note: If MXLEAD > 1, forecasts are returned for t = NOBS + 1, .NOBS + 2, …, NOBS + MXLEAD
Default: MXLEAD = 0 (no forecasts).

MAXIT — Maximum number of iterations. (Input)
Default: MAXIT = 50.

IPRINT — Printing option. (Input)

 

IPRINT	Action
0	No printing
1	Prints final results only
2	Prints intermediate and final results

Default: IPRINT = 0.

PREG — Array of length K + t containing the estimated regression coefficients, where t=0 if TREND=.FALSE. or t=1 if TREND=.TRUE.. (Output)

SE — Array of length p + q containing the standard errors of the ARMA parameter estimates. (Output)

AVAR — White noise variance estimate. (Output)
Note: If IMODEL(0)+IMODEL(2)= 0 and K>0, AVAR is the mean squared regression residual.

REGSE — Array of length K + t containing the standard errors of the regression estimates, where t=0 if TREND=.FALSE. or t=1 if TREND=.TRUE.. (Output)

PREGVAR — Array of length (K + t) by (K + t) containing the variances and covariances of the regression coefficients, where t=0 if TREND=.FALSE. or t=1 if TREND=.TRUE.. (Output)

AIC — Akaike's Information Criterion for the fitted ARMA model. (Output)

LLIKE — Value of –2(ln(likelihood)) for fitted model. (Output)

FCST — Array of length MXLEAD containing the forecasts for time points t = NOBS + 1, .NOBS + 2, …, NOBS + MXLEAD. (Output)

FCSTVAR — Array of length MXLEAD containing the forecast variances for time points t = NOBS + 1, .NOBS + 2, …, NOBS + MXLEAD. (Output)

Fortran 90 Interface

Generic: CALL REG_ARIMA (Y, IMODEL, PARMA [, …])

Specific: The specific interface names are S_REG_ARIMA and D_REG_ARIMA.

Description

Routine REG_ARIMA fits an ARIMA(p, d, q) to a univariate time series with the possible inclusion of one or more regression variables.

Suppose Yt, t = 1, …, N, is a time series such that the d‑th difference is stationary. Further, suppose at is a series of uncorrelated, mean 0 random variables with variance .

The Auto‑Regressive Integrated Moving Average (ARIMA) model for {Yt, at} can be expressed as

 

where B is the backshift operator,

 

and

 

The notation for this model is ARIMA(p, d, q) where p is the order of the autoregressive polynomial , d is the order of the differencing needed to make stationary, and q is the order of the moving‑average polynomial .

The ARIMA model can be extended to include regression variables , by using the residuals (of the multiple regression of on ) in place of in the above ARIMA model.

 

Equivalently,

 

where

 

is the differenced residual series.

To estimate the (p + q + K) parameters of the specified regression ARIMA model, REG_ARIMA uses the iterative generalized least squares method (IGLS) as described in Otto, Bell, and Burman (1987).

The IGLS method iterates between two steps, one step to estimate the regression parameters via generalized least squares (GLS) and the second step to estimate the ARMA parameters. In particular, at iteration m, the first step finds

 

by solving the GLS problem with weight matrix

 

where

 

That is, minimizes , where , is an N by K matrix with i‑th column, , and is an N by N weight matrix defined using the theoretical autocovariances of the series . The series is modeled as an ARMA(p,q) with parameters and . At iteration m, the second step is then to obtain new estimates of and for the updated series,. To find the estimates and , REG_ARIMA uses the exact likelihood method as described in Akaike, Kitagawa, Arahata and Tada (1979) and used in routine, MAX_ARMA.

Comments

When forecasts are requested (MXLEAD > 0), REG_ARIMA requires that future values of the independent variables are provided in optional argument XLEAD. In effect, REG_ARIMA assumes the future X’s are known without error, which is valid for any deterministic function of time such as a seasonal indicator. Also, in economics, certain factors that are considered to be leading indicators are treated as deterministic for the purpose of predicting changes in the economy. Users may consider using a more general transfer function model if this is an unreasonable assumption. REG_ARIMA calculates forecast variances using the asymptotic result found in Fuller (1996) , Theorem 2.9.4. To obtain the standard errors of the ARMA parameters, REG_ARIMA calls routine NSLSE for the final w series.

Examples

Example 1

The data set consists of annual mileage per passenger vehicle and annual US population (in 1000’s) spanning the years 1980 to 2006 (U.S. Energy Information Administration, 2008). Consider modeling the annual mileage using US population as a regression variable.

 

use reg_arima_int

use umach_int

implicit none

 

integer, parameter :: mxlead=5, idep=2

integer :: i, nout, nobs, n

integer :: imodel(3)=(/1,0,0/), indind(1)=(/1/)

real(kind(1e0)) :: y(29), parma(2), xlead(mxlead, 2)

real(kind(1e0)) :: preg(2), regses(2)

real(kind(1e0)) :: ses(1), fcst(mxlead), fcstvar(mxlead)

real(kind(1e0)) :: avar,llike

! US mileage and population (1980-2008)

! Source: U.S. Energy Information

! Administration (Oct 2008).

 

real(kind(1e0)) :: x(29,2)

data x / &

22722.4681, 22946.5714, 23166.4458, 23379.1990, &

23582.4902, 23792.3795, 24013.2887, 24228.8918, &

24449.8982, 24681.923, 24962.2814, 25298.0941, 25651.4224,&

25991.8588, 26312.5820999999, 26627.8393, 26939.4284, &

27264.6925, 27585.4104, 27904.0168, 28217.1936, &

28503.9803, 28772.6647, 29021.0914, 29289.2127, &

29556.0549, 29836.2973, 30129.0332, 30405.9724, &

9062.0, 8813.0, 8873.0, 9050.0, 9118.0, 9248.0, 9419.0, &

9464.0, 9720.0, 9972.0, 10157.0, 10504.0, 10571.0, &

10857.0, 10804.0, 10992.0, 11203.0, 11330.0, 11581.0, &

11754.0, 11848.0, 11976.0, 11831.0, 12202.0, 12325.0, &

12460.0, 12510.0, 12485.0, 12293.0/

 

call umach(2,nout)

! Example 1

! The first column is the scaled US

! population and the second column is

! the annual mileage per vehicle

n = size(x,1)

nobs = n - mxlead

call scopy(nobs,x(1:n,idep),1,y,1)

call reg_arima(y, imodel ,parma, nobs=nobs, iprint=1, x=x, &

indx=indind, avar=avar, llike=llike, preg=preg, &

regse=regses, mxlead=mxlead, xlead=x(nobs+1:n,1:2), &

trend=.true., fcst=fcst, fcstvar=fcstvar, se=ses)

end

Output

 

Final results for regarima model (p,d,q) =  1 0 0

 

 Final AR parameter estimates/ std errors

 

          0.73063        0.13499

 

-2*ln(maximum log likelihood) =   231.8354

 

 White noise variance 10982.5654

 

 Regression estimates:

 

               coef          reg. se

     1      -3481.22607      689.02661

     2          0.54237        0.02673

 

 Forecasts with standard deviation

 

 t        Y fcst          Y fcst std dev

 25      12360.40137        153.89659

 26      12514.61035        171.89198

 27      12673.53320        180.76634

 28      12837.36719        185.32974

 29      12991.26855        187.72037

Example 2

The data set consists of simulated weekly observations containing a strong annual seasonality. The seasonal variables are constructed and sent into REG_ARIMA as regression variables.

 

use reg_arima_int

use umach_int

use const_int

implicit none

 

integer, parameter :: mxlead = 4

integer :: nobs, i,n,idep,nout

integer :: imodel(3) =(/2,0,0/)

real(kind(1e0)) :: PI

real(kind(1e0)) :: preg(3),regses(3),parma(3)

real(kind(1e0)) :: ses(2),fcst(mxlead),fcstvar(mxlead)

real(kind(1e0)) :: avar,llike,aic,tmpr

real(kind(1e0)) :: x(100,2), xlead(mxlead,2)

real(kind(1e0)) :: y(104) =(/ &

32.27778,32.63300,33.13768,34.4517,34.63824, &

37.31262,37.35704,37.03092,36.39894,35.75541,&

35.10829,34.70107,34.69592,32.75326,30.85370,&

31.10936,29.47493,29.14361,28.50466,30.09714,&

28.49403,27.23268,23.49674,22.71225,21.42798,&

18.68601,17.40035,16.06832,15.31862,14.75179,&

13.40089,13.01101,12.44863,11.27890,11.51770,&

14.31982,14.67036,14.76331,15.35644,17.04353,&

18.39931,18.21919,18.72777,19.61794,22.31733,&

23.79600,25.41326,25.60497,27.93579,29.21765,&

29.60981,28.46994,28.78081,30.96402,35.49537,&

35.75124,36.18933,37.2627,35.02454,33.57089, &

35.00683,34.83886,34.19827,33.73966,34.49709,&

34.07127,32.74709,31.97856,31.3029,30.21916, &

27.46015,26.78431,25.32815,23.97863,21.83837,&

21.00647,20.58846,19.94578,17.38271,17.12572,&

16.71847,17.45425,16.15050,13.07448,12.54188,&

12.42137,13.51771,14.84232,14.28870,13.39561,&

15.48938,16.47175,17.62758,16.57677,18.20737,&

20.8491,20.15616,20.93857,23.73973,25.30449, &

26.51106,29.43261,32.02672,32.18846/)

 

 

call umach(2,nout)

PI = const('PI')

 

! The data are simulated weekly

! observations with an annual

! seasonal cycle

n = size(y)

nobs = n - mxlead

! Create the seasonal variables

do i = 1, nobs

x(i,1)=sin(2*PI*i/float(52))

x(i,2) = cos(2*PI*i/float(52))

end do

do i=1, mxlead

xlead(i,1)=sin(2*PI*(i+nobs)/float(52))

xlead(i,2) = cos(2*PI*(i+nobs)/float(52))

end do

 

call reg_arima(y, imodel, parma, iprint=1, x=x, &

nobs=nobs, avar=avar, llike=llike, &

preg=preg, regse=regses, mxlead=mxlead, &

xlead=xlead, trend=.true., fcst=fcst, &

fcstvar=fcstvar, se=ses)

end

Output

 

Final results for regarima model (p,d,q) =  2 0 0

 

 Final AR parameter estimates/ std errors

 

          0.71860        0.09836

 

         -0.25991        0.09827

 

-2*ln(maximum log likelihood) =   -13.6209

 

 White noise variance     0.8783

 

 Regression estimates:

 

               coef          reg. se

     1         24.81010        0.17175

     2          9.68013        0.23993

     3          5.72305        0.24751

 

 Forecasts with standard deviation

 

 t        Y fcst          Y fcst std dev

101         26.74492          1.31358

102         28.07805          1.47616

103         29.33707          1.49560

104         30.53160          1.49560