|
Fits a univariate, non-seasonal ARIMA time series model with the inclusion of one or more regression variables.
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)
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, in
which case
K = size(INDX) .
Default: No
regression variables are included unless TREND = .TRUE., then a trend
variable is included.
XLEAD Array of
size MXLEAD by
K containing the regression data to be used in obtaining forecasts, where
K = size(X,2) (or K=
size(INDX) if INDX is supplied) is
the number of user supplied regression variables.
(Input)
Note: If MXLEAD > 0 and
optional argument X is present, XLEAD is required.
INDX Index array containing the column numbers in X that are to be used for the regression variables. (Input)
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 .
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 . (Output)
FCSTVAR Array
of length MXLEAD
containing the forecast variances for time points . (Output)
Generic: CALL REG_ARIMA (Y, IMODEL, PARMA [, ])
Specific: The specific interface names are S_REG_ARIMA and D_REG_ARIMA.
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 ,
, is a time series such that the d-th difference is
stationary. Further, suppose
is a series of
uncorrelated, mean 0 random variables with variance
.
The
Auto-Regressive Integrated Moving Average (ARIMA) model for 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.
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 Xs 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.
The data set consists of annual mileage per passenger vehicle and annual US population (in 1000s) 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
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
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
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
PHONE: 713.784.3131 FAX:713.781.9260 |