Short Term Speed Variance Forecasting Using Linear Stochastic Modeling of Univariate Traffic Speed Series
Jianhua Guo
Research Associate
Brian L. Smith
Associate Professor
University of Virginia
Center for Transportation Studies
July 16, 2007
A Research Project
Report for the Intelligent Transportation Systems Implementation Center (ITS)
A U.S. DOT University Transportation Center
Dr. Brian L. Smith, PhD
Department of Civil Engineering
Email: bls2z@virginia.edu
Center for Transportation Studies at the University of Virginia produces outstanding transportation professionals, innovative research results and provides important public service. The Center for Transportation Studies is committed to academic excellence, multi-disciplinary research and to developing state-of-the-art facilities. Through a partnership with the Virginia Department of Transportation’s (VDOT) Research Council (VTRC), CTS faculty hold joint appointments, VTRC research scientists teach specialized courses, and graduate student work is supported through a Graduate Research Assistantship Program. CTS receives substantial financial support from two federal University Transportation Center Grants: the Mid-Atlantic Universities Transportation Center (MAUTC), and through the National ITS Implementation Research Center (ITS Center). Other related research activities of the faculty include funding through FHWA, NSF, US Department of Transportation, VDOT, other governmental agencies and private companies.
Disclaimer: The contents of this report reflect the views of the authors, who are responsible for the facts and the accuracy of the information presented herein. This document is disseminated under the sponsorship of the Department of Transportation, University Transportation Centers Program, in the interest of information exchange. The U.S. Government assumes no liability for the contents or use thereof.
|
1. Report No. |
2. Government Accession No. |
3. Recipient’s Catalog No. |
|||||
|
UVACTS-15-17-10 |
|
|
|||||
|
4. Title and Subtitle |
5. Report Date |
||||||
|
Short Term Speed Variance Forecasting Using Linear
Stochastic Modeling of Univariate Traffic Speed Series |
July 17, 2007 |
||||||
|
|
6. Performing Organization Code |
||||||
|
|
|
||||||
|
7. Author(s) |
8. Performing Organization Report No. |
||||||
|
Dr. Brian Smith, Jianhua Guo |
|
||||||
|
9. Performing Organization and Address |
10. Work Unit No. (TRAIS) |
||||||
|
Center for Transportation Studies |
|
||||||
|
University of Virginia |
11. Contract or Grant No. |
||||||
|
PO Box 400742 Charlottesville, VA 22904-7472 |
|
||||||
|
12. Sponsoring Agencies' Name and Address |
13. Type of Report and Period Covered |
||||||
|
Office of University
Programs, Research Innovation and Technology Administration US Department of
Transportation 400 Seventh Street, SW Washington DC 20590-0001 |
|
Final Report |
|||||
|
|
|
14. Sponsoring Agency Code |
|||||
|
|
|
|
|||||
|
15. Supplementary
Notes |
|||||||
|
|
|||||||
|
16. Abstract |
|||||||
|
|
|||||||
|
17 Key Words |
18. Distribution Statement |
||||||
|
Speed Variance Forecasting, Univariate Traffic
Speed Series |
No restrictions. This document is available to the
public. |
||||||
|
19. Security Classif. (of this report) |
20. Security Classif. (of this page) |
21. No. of Pages |
22. Price |
||||
|
Unclassified |
Unclassified |
21 |
N/A |
||||
1 Introduction
Intelligent transportation systems (ITS) offers the potential to address critical transportation needs. However, most ITS currently operates in a reactive mode. While this provides some level of benefit, “the full benefits of ITS cannot be realized without an ability to anticipate traffic conditions in the short-term (less than one hour into future).” (Smith and Oswald, 2003) Based on the anticipated traffic condition, proactive transportation management and comprehensive traveler information services are feasible. Therefore, traffic condition forecasting has been identified as one of the major challenges for ITS.
Considering the forecasting process as a state extrapolation process governed by certain regularity, the development of traffic condition forecasting methods demands a sound understanding of traffic condition dynamics. Volume, speed, and density (or occupancy for the widely-deployed inductive loop detectors) are three traffic variables that are most commonly used to characterize traffic conditions, and suitable traffic condition forecasting methods are expected to be built upon traffic condition dynamics in terms of these traffic variables.
Previous efforts addressing traffic flow forecasting at a higher aggregation interval, such as 15-minutes, indicate traffic conditions to be linear stochastic. The traffic flow forecasting methods can be roughly classified into nonlinear theory based methods and linear theory based methods. The former assumes traffic dynamics are nonlinear, and can be emulated through nonlinear operations. Typically, this category includes non-parametric regression, neural networks, kernel smoothing, and local linear regression. The latter assumes traffic dynamics are linear and can be emulated through linear operations. Typically, this category includes univariate Box-Jenkins approach, exponential smoothing, spectral analysis, and multivariate time series methods. Adaptive methods, such as Kalman filter and recursive least square, can be classified into linear methods due to their nature as sequential projection in linear space. Williams (1999), Smith et al. (2002), and Williams and Hoel (2003) showed that traffic flow forecasting method based on Seasonal Autoregressive Integrated Moving Average (SARIMA) process outperformed nonlinear theory based methods, supporting the adoption of SARIMA process to describe traffic flow dynamics. Guo (2005) further appended a Generalized Autoregressive Conditional Heteroscedasticity (GARCH) model to describe the conditional variance for forecasting confidence interval construction. In addition, recent research on traffic state transition summarized in Daganzo (2002) revealed that traffic dynamics could be modeled using simple first order continuum theory, suggesting a linear regularity in traffic dynamics.
Based on linear traffic dynamics, the traffic speed series is naturally expected to be described using linear stochastic theory. In this report, first, the conditional mean of traffic speed series will be investigated using ARIMA model. Second, the conditional variance of traffic speed series will be investigated using GARCH model. On one hand, forecast confidence intervals could be constructed using conditional variance computed from GARCH process, indicating the certainty associated with the forecasts (or predictability of traffic). On the other hand, conditional variance could be utilized in more broad ITS settings, one of which is to determine the sample size for probe vehicle based traffic surveillance systems. In summary, the combination of ARIMA model and GARCH model will present a complete structure describing the first two conditional moments of traffic speed series, facilitating the development of an online forecasting system.
The rest of this report is organized as follows. Section 2 discusses the related literature; section 3 presents an overview of time series modeling; section 4 applies the model to the real traffic speed series from the US and the UK; finally, section 5 concludes this report with a discussion of the results.
2 Related Literature
Like the efforts on traffic flow forecasting that indicate traffic flow series can be best described as SARIMA process (Williams 1999; Smith et al., 2002; Williams and Hoel, 2003), the investigation of traffic speed series was largely contained within the development of traffic speed estimation and forecasting methods.
§ Yang et al. (2004) provided an online recursive traffic speed forecasting method based on an arbitrary autoregressive process. The authors use Kalman filter to process this autoregressive model and provide forecasting confidence interval using prediction variance, which will converge to a constant eventually. Actually, the maximum difference between the lengths of confidence intervals (Figure 6 in their paper) is about 0.1 mph for confidence interval width around 8 mph.
§ Based on nonlinear dynamics assumption, Sun et al. (2003) used local linear regression method to forecast 5-minute traffic speed data collected on the US-290 Northwest freeway in Houston. A comparative study showed that local linear regression outperformed k-nearest neighbor and kernel smoothing method. In addition, Sun et al. (2004) presented a method to construct prediction confidence intervals for predictions from local linear regression. Asymptotic prediction confidence intervals and prediction confidence intervals based on bootstrap method were constructed and the bootstrap method was recommended. Their results showed that a large proportion of real observations were outside of prediction confidence intervals given by the bootstrap method, which invalidated the goodness of their approach.
§
Different from univariate analysis, Tavana et al.
(2000) approached using transfer function model (bivariate time series model).
In their work, based on the concept of dynamic speed-density relationship, the
variation of speed around a static speed-density relation is captured by a
transfer function model, in which density (estimated from occupancy) was used
as forcing function and the deviation of speed from its static equilibrium
value was used as system output. The static equilibrium speed-density relation
was assumed to be linear and was estimated from speed-density data. The speed
equilibrium value was then computed using density and the estimated static
equilibrium speed-density relation. First-order differencing was performed on
the input and output to introduce stationarity. The transfer model was
identified using cross-correlation analysis after applying pre-whitening
operation to both input and output series. The noise term after removing the
effects of cross-correlation was modeled as MA(1). In essence, the transfer
function model is a regression of the output (speed deviation) onto input
(density) with MA(1) errors. The model parameters were estimated through an
optimization process minimizing the conditional sum square of errors. In
forecasting speed for time 
, the density at time 
needs to be obtained
first. Two options were suggested: prediction using univariate time series
forecasting method and simulation. The latter was adopted in their work.
§ Based on the transfer function model in Tavana et al. (2000), Huynh et al. (2002) proposed an adaptive process of estimating the model parameters. Results from simulation-based analysis showed that adaptive model outperformed non-adaptive model. Qin et al. (2004) proposed a calibration of the same model using real sensor data from Irvine, CA, and reaffirmed that adaptive model outperformed the non-adaptive model and also showed that the increase of number of links applied in the model improved speed estimation.
§ Antoniou et al. (2005) approached using state space model. The speed-density relation was regarded as the measurement equation, and the parameters in the equation were regarded as system state. The transition of system state was formulated as an autoregressive process. The degree of the transition process, the transition matrix, and the errors were estimated off-line. Three variations of Kalman filter algorithm, including extended Kalman filter (EKF), Iterated EKF, and Unscented Kalman filter (UKF), were tested using real sensor data from Europe and the US. Numerical results showed that speed estimation using online calibrated speed-density relation was improved over speed estimated from off-line calibrated speed-density relation. As is in the transfer function model (Tavana et al. 2000, Huynh et al. 2002, and Qin et al. 2004), the prediction of the occupancy (density) has to be obtained before obtaining the forecast of speed.
3 Time Series Models
In this section, the theory of
time series modeling is briefly reviewed for completeness (Fuller 1996). A time series defined as 
can be strictly
stationary or weakly stationary. Given 
as the probability
distribution function and lag 
, a time series is strictly stationary if
(1)
for all the
possible sets of 
and
. This definition states that the joint probability of any
given subset of a strictly stationary time series is independent of the lag. A time
series is weakly stationary if
(1) the expected value of 
is constant for all 
; and
(2) the covariance matrix of 
is the same as that
of 
for all finite sets 
and 
.
Weakly stationarity only
restricts the first and second moment of a time series, and hence it is also
termed as covariance stationarity. According to the Wold Decomposition Theorem,
a weakly stationary time series 
defined on 
can be represented as
,
(2)
where 
is a deterministic
component and 
can be modeled as an
autoregressive moving average (ARMA) process defined as
,
(3)
where
: backshift operator such that 
;
: short-term AR polynomial with order 
;
: short-term MA polynomial with order 
;
For a non-stationary time series
, the non-stationarity can be handled by proper differencing,
and the whole model, seasonal integrated autoregressive moving average-SARIMA (p,d,q)(P,D,Q)S , is defined as
(4)
where
: seasonal differencing with seasonality 
and order 
;
: normal differencing
with order 
;
: 
order seasonal AR polynomial;
: 
order seasonal MA polynomial.
It is assumed that the roots of 
, 
, 
, and 
are outside the unit
circle, which enables the causality and invertibility of the process; also, these
polynomials should have no common
factors. 
is assumed to be Gaussian
white noise with mean zero and constant variance 
; in addition, it is further assumed to follow GARCH process
defined as below.
, 
, 
(5)
where
: conditional
variance at time 
, i.e., 
with 
as the information up
to time 
;
,
: order of the
GARCH process with 
and 
;
: positive
constant coefficient;
: non-negative
coefficients of the lagged sample variance 
;
: non-negative
coefficients of the lagged conditional variance 
;
If 
, GARCH(p,q) process is reduced to ARCH (
) process; if 
, the GARCH(p,q) process is simplified into Gaussian
white noise, and the SARIMA+GARCH structure is reduced to SARIMA process. An
alternative representation shows that GARCH process can be interpreted as ARMA
process in 
as:
(6)
where 
and 
(Bollerslev 1986), indicating
that the modeling of GARCH model is similar to that of SARIMA in the sense of
squared series. Box-Jenkins approach provided a complete framework of
SARIMA modeling, including identification, model estimation, forecasting, and
diagnostic check. Interested readers are referred to (Box and Jenkins 1994) for
detail.
The combination of equation (4) and (5) presented a complete stochastic
structure, SARIMA+GARCH. Using results in Engle (1982) and Bollerslev (1986),
this structure can be processed separately without efficiency lost: (1) the
processing of SARIMA, and (2) the processing of GARCH using the residuals from
the SARIMA model.
4 Empirical Study
This section applied the SARIMA+GARCH structure to traffic speed series. First, the sensor data used in this analysis is discussed; second, using the sensor data, the model for speed levels is identified; third, heteroscedasticity tests and GARCH model selection are discussed; finally, the adequacy of the identified stochastic structure is demonstrated by investigating its performance.
4.1 Data
Taffic data collected from two highway systems, including the motorway system in the United Kingdom and freeway systems in Richmond, Virginia in the US, are used in this study. In this study, these traffic speed data were aggregated into 15-minute data according to Edie (1963), i.e., flow weighted harmonic mean. Obvious erroneous data are excluded and treated as missing values. In total, 10 traffic speed series were compiled (see Table 1).
Table
1: Overview of Selected Traffic Speed
Series
|
Region |
Station |
Total Length |
Percent
Missing(%) |
|
RHMD |
12040 |
6720 |
7.05 |
|
RHMD |
12042 |
6666 |
5.55 |
|
RHMD |
12043 |
6666 |
4.64 |
|
UK |
2737a |
30912 |
4.71 |
|
UK |
2808b |
20676 |
5.05 |
|
UK |
4565a |
6289 |
4.55 |
|
UK |
4762a |
35039 |
3.93 |
|
UK |
4897a |
30912 |
7.87 |
|
UK |
6951a |
35039 |
1.78 |
|
UK |
6954b |
35039 |
1.92 |
Note: UK: the United Kingdoms; RHMD: Richmond, VA.
4.2 ARIMA Model Selection
In this section, the order of ARIMA model is determined. Unlike traffic flow series in which strong daily or weekly pattern exists (Williams 1999), for speed series, the daily/weekly seasonality is not significant. According to previous research, the speed remains nearly constant even at quite high flow rates (Revised Traffic Flow Theory Monograph, 1994). The same argument is also reflected in Highway Capacity Manual 2000 (See Exhibit 8-26 and Exhibit 8-27 for speed and volume variations across weekdays and Saturday, and Exhibit 13-2 for the speed-flow relationships for the basic freeway segment). All these evidences lead to the fact that the seasonal factor of the SARIMA process could be eliminated in modeling traffic speed series.
Although speed series will roughly maintain at a level for most of the time, speeds do oscillate due to disturbances like incident and weather, causing non-stationarity in the series. In this study, the speed series will be differenced once to introduce stationarity, which is in accordance with previous efforts in (Tavana et al. 2000, Huynh et al. 2002, and Qin et al. 2004).
After the first order differencing, and ignoring the seasonality, an autoregressive moving average (ARMA) model could be constructed. The sample ACF and sample PACF from lag 1 to lag 16 for the differenced series are computed using SAS ®. See Figure 1 below.
RHMD
12040 PACF RHMD
12040 ACF
RHMD
12042 PACF RHMD
12042 ACF
RHMD
12043 PACF RHMD
12043 ACF
UK
2737a PACF UK
2737a ACF
UK
2808b PACF UK
2808b ACF
UK
4565a PACF UK
4565a ACF
UK
4762a PACF UK
4762a ACF
UK
4897a ACF UK
4897a PACF
UK
6951a PACF UK
6951a ACF
UK
6954b PACF UK
6954b ACF
Figure 1:
Sample ACF and Sample PACF of First-Order Differenced Speed Series
(x-axis: lag; y-axis: ACF or PACF; Dashed line around x-axis: 95% confidence
limits)
The ACF (cutoffs at lag one) and PACF (tails off) for Richmond stations indicates the series can be modeled as MA(1) process. The ACF and PACF for UK stations are quite volatile, and no single AR model or MA model could be used. However, one common point is that the magnitude of the ACF and PACF of these stations is small. Under this circumstance, combing the consideration of Richmond station, in this study, we use ARMA(1,1) to model the first-order differenced speed series for all the stations. To summarize, the speed level process is modeled as ARIMA(1,1,1). It is worthwhile to note that the model is meant to describe the general trend of speed process in the long run for building an online forecasting model, and this model will not necessarily be the exact model for each series. The discrepancy will be compensated in the online forecasting algorithm that can adapt to local variations.
4.3 Heteroscedasticity Test and GARCH Model Selection
Upon the Selection of ARIMA(1,1,1), the modeling of the second conditional moments is described as follows. First, missing values were replaced by the one-step-ahead forecasts using ARIMA(1,1,1) model in an iterative procedure. The iteration will stop when the model standard error stops decreasing or the maximum number of iterations has been reached. Based on the missing value imputed series, ARIMA(1,1,1) model was estimated and the forecasting residuals were computed. In this process, the model parameters were estimated using maximum likelihood method. Then, using the forecasting residuals, heteroscedasticity test using Langrage multiplier test and portmanteau test was performed to justify the addition of the GARCH process. To account for test inflation due to series length, four strategies, i.e., test by length, by month, by week, and by day, are applied. In test by length, the entire residual series is tested; for the other three tests, the residual series is first broken into sub-series and then tested individually. According to the length of series in each test, the testing power is decreasing in the order of test by length, month, week, and day. To be conservative, for each series or sub-series, statistical significant heteroscedasticity will be claimed only when both tests (Langrage multiplier test and portmanteau test) are significant at <0.001 level for lags from 1 to 12. For test by month, week, and day, the percentages of the significant sub-series will be computed to show the degree of variance change for each station.
The results for test by month/week/day are presented in Table 2. For test by length, the results are uniformly significant across stations. It can be seen that the process variance is changing across time, and the percentage of significance for each test decreases with the decreasing of testing power.
Table 2: Heteroscedasticity Test Results
|
Region |
Station |
|
Test by Month |
|
Test by Week |
|
Test by Day |
||||||
|
|
A |
B |
% |
|
A |
B |
% |
|
A |
B |
% |
||
|
RHMD |
12040 |
|
3 |
2 |
67 |
|
10 |
2 |
20 |
|
70 |
2 |
3 |
|
RHMD |
12042 |
|
3 |
3 |
100 |
|
10 |
2 |
20 |
|
70 |
0 |
0 |
|
RHMD |
12043 |
|
3 |
3 |
100 |
|
10 |
10 |
100 |
|
70 |
2 |
3 |
|
UK |
2737a |
|
11 |
10 |
91 |
|
46 |
19 |
41 |
|
322 |
10 |
3 |
|
UK |
2808b |
|
8 |
7 |
88 |
|
31 |
22 |
71 |
|
216 |
7 |
3 |
|
UK |
4565a |
|
3 |
3 |
100 |
|
10 |
9 |
90 |
|
66 |
3 |
5 |
|
UK |
4762a |
|
12 |
12 |
100 |
|
53 |
45 |
85 |
|
365 |
56 |
15 |
|
UK |
4897a |
|
11 |
11 |
100 |
|
46 |
36 |
78 |
|
322 |
17 |
5 |
|
UK |
6951a |
|
12 |
11 |
92 |
|
53 |
41 |
77 |
|
365 |
13 |
4 |
|
UK |
6954b |
|
12 |
12 |
100 |
|
53 |
42 |
79 |
|
365 |
20 |
5 |
Note: A – total number of month/week/day; B – number of significant month/week/day.
According to equation (6), GARCH process can be regarded as ARMA process in terms of squared series; therefore, the identification can be performed through investigating the sample ACF and PACF of squared residual series. In this study, GARCH(1,1) is selected, and as the case for ARMA(1,1) for first-order differenced speed series, GARCH(1,1) is used to describe the general trend of the changing variance across stations and the discrepancy can be compensated for when developing the online algorithm.
4.4 Model Adequacy Evaluation
Summarizing section 4.2 and 4.3, the model for traffic speed series was identified as ARIMA(1,1,1)+GARCH(1,1), which is presented below.
(7)
, 
, 
(8)
In this subsection, the adequacy of this structure will be demonstrated through investigating its performance.
4.4.1 Performance Measures
The
performance measures include the measures of forecasting accuracy and
prediction confidence intervals. Given 
as observations, 
as forecasts, and 
as the total number
of observations, the forecasting accuracy measures are defined as: (1) Mean
Absolute Error: 
; (2) Mean Absolute Percentage Error: 
; (3) Root Mean
Square Error: 
. The evaluation of
confidence intervals is complicated with multiple criteria. In this study,
following two intuitions, i.e., (1) real observations are expected to fall
within the prediction confidence intervals, and (2) narrower prediction
confidence intervals are more informative than wider ones, two summary
statistics, e.g., prediction
confidence interval kickoff percentage and average prediction confidence
interval width to speed ratio are constructed. A prediction confidence interval kickoff is identified as an
observation falling outside of the corresponding prediction confidence
interval, and a percentage can be calculated as the total number of kickoffs
divided by the total number of observations. Discussed in (Kang and Schmeiser
1990; Schmeiser and Yeh 2002), suitable prediction confidence intervals will
have a kickoff percentage close to the nominal value, i.e., 5% for 95%
confidence interval, and under this condition, the narrower the interval, the
better. In computing performance
measures, only non-missing speeds and corresponding forecasts and prediction
confidence intervals are used. The imputed speeds are excluded from the
computation.
4.4.2 Forecasting Accuracy
The forecasting accuracy measures are listed in Table 3. It can be seen the forecasting accuracy measures of all the stations (except for station RHMD 12043) are satisfactory: the forecasting error generally varies around 1.5-2 mph, with relative error around 2.5-4%. The high forecasting performance is in accordance with the speed series dynamics. As discussed previously, the speed will generally maintain at a fixed level except for occasional period when traffic is subjected to external disturbances. For RHMD 12043, the forecasting measures indicate that speed series for this station is not as predictable as other stations. This argument will be further substantiated when investigating the forecasting confidence interval.
Table 3: Forecasting Accuracy Measures
|
Region |
Station |
MAE |
MAPE |
RMSE |
|
RHMD |
12040 |
2.17 |
4.07 |
3.59 |
|
RHMD |
12042 |
1.63 |
3.00 |
2.55 |
|
RHMD |
12043 |
7.56 |
22.97 |
10.58 |
|
UK |
2737a |
1.38 |
2.23 |
2.20 |
|
UK |
2808b |
1.65 |
3.04 |
3.04 |
|
UK |
4565a |
1.50 |
2.96 |
2.89 |
|
UK |
4762a |
2.32 |
4.85 |
4.36 |
|
UK |
4897a |
1.55 |
2.75 |
2.69 |
|
UK |
6951a |
1.46 |
2.75 |
2.52 |
|
UK |
6954b |
1.40 |
2.71 |
2.87 |
4.4.3 Prediction Confidence Interval
The performance measures of prediction confidence intervals are listed in Table 4. It can be seen that the prediction confidence interval kickoff percentage for all the stations varies around 3-5%, indicating the proposed prediction confidence interval construction method based on GARCH model is workable although the results is slightly conservative in yielding an high coverage. In the mean time, the average prediction confidence interval width to speed ratio are around 0.15-0.25, which can be translated to 9-15 mph using an assumed speed level at 60 mph. This is approximately ±5 –±8 mph around traffic speed level. The average prediction confidence interval width to speed ratio for RHMD 12043 is 0.93, indicating a wide confidence interval for the forecasts for this station and thus a low predictability. This is in accordance with the observation on the forecasting accuracy measures for this same station in Table 3.
