Estimation and Accuracy of Origin-Destination Highway Freight Weight and Value Flows

Paul Metaxatos
University of Illinois at Chicago

ABSTRACT

This paper proposes a spatial interaction modeling framework and implements a maximum likelihood estimation of highway freight weight and value flows using the gravity model. The computation of the standard error of the flow estimates provides the basis for measuring the level of accuracy of the estimates. The results provide evidence of the suitability of gravity models for freight forecasting given the excellent fit and the small variances.

INTRODUCTION

The measurement of freight movements requires tracking freight flows across geographic and political boundaries. This is a particularly challenging task given the current capabilities for state and regional data acquisition. Various mathematical approaches have been implemented (Memmot 1983; USDOT 1996; Cambridge Systematics 1997) to circumvent this problem, but none, to the best of this author's knowledge, proposes a measure to assess the accuracy of the computed flows.

This paper proposes to fill this gap using developments in spatial interaction modeling that have not been demonstrated on a large scale to date. The methodology computes maximum likelihood flow estimates and obtains their covariance matrices that, in turn, may be used to obtain confidence intervals and carry out certain tests of hypotheses. The approach can accommodate the large number of origins and destinations typically encountered in freight (and passenger) travel forecasting.

The methodology was applied to highway freight weight and value flows of international trade traffic between seaports or border ports and destination states (see Metaxatos (2002) for details). The variances of the flow estimates computed were remarkably small. The demand for freight transportation flows can then be estimated within a desired confidence level. Moreover, the empirical analysis undertaken provides evidence that the theoretical framework proposed in this paper is rich enough for freight demand forecasting applications.

THEORETICAL FRAMEWORK

Commodity shipments in this paper are thought to be realized patterns of spatial interactions that typically result from many independent decisions by individual firms, each constituting a relevant subsystem within the economy as a whole. Hence, if the travel behavior of each firm is modeled as a very small interaction process, the resultant interaction process can be taken to be the superposition of all these processes. It may be argued that for large collections of small frequency processes, the resulting superimposed process is approximately Poisson and, therefore, completely characterized by its associated mean interaction frequencies (Sen and Smith 1995).

In this light, assuming that the observations N_ij of shipment weight and value between origin seaports/border ports of entry i and destination states j can be described by the gravity model, then

(1)

In this paper, T_ij 's (the stochastic term) are interpreted as the expected international trade traffic flow (in terms of weight and value) carried by highway from external station i to state j. The A_i 's are factors related to the origin zone i and the B _j ' s are destination-related factors. The F_ij 's are factors that reflect the separation between i and j. A common form that is general enough for most applications is

(2)

This form is called an exponential form and are different measures of separation, while θ_k 's are parameters to be estimated. Potential measures of separation include travel time, distance, generalized costs, etc.

In the gravity model, observable quantities N_i* = Σ_j N_ij,N_*j = Σ_i N_ij,N_ij and their expected values T_i*, T_*j, and T_ij are described by means of an underlying structure consisting of unobservable quantities A_i, B_j, and F_ij. Similar situations abound in statistics. In moving average models, for example, observations are described by means of unobservable parameters. A like situation exists in analysis of variance models.

Although the origin and destination factors are unobservable, they do have physical interpretations. For example, if for some origin i, there are two destinations j and j ' such that F_ij = F_ij', then T_ij /T_ij' = B_j /B_{j '}. Thus other factors being equal, T_ij is proportional to B_j (but, in general, not proportional to T_*j), and is called the attractiveness of j. Similarly, the origin factor A_i may also be called the emissiveness of i.

Clearly, BjF_ij is the effect of the destination factor Bj at i , or the accessibility of j as perceived from i . This is a spatial analogy of the temporal concept of present value in economics, where a dollar earned in n years in the future is worth only (1 + σ)^–n now, where σ is the interest rate. Similarly, A_iF_ij is the effect of the origin factor A_i at j . The sum α_i = Σ_jB_jF_ij may be called the total accessibility of all destinations at i , and the sum β_j = Σ_iA_iF_ij may be called the total accessibility of all origins at j . If, for example, T_i* is kept fixed as α_i increases, the push A_i decreases. Thus, as the competition α_i from the destinations increases, the push at i decreases. From the point of view of someone at i , α_i measures accessibility; from the viewpoint at j , it measures competition. Similar statements can be made about B_j.

Maximum Likelihood Estimation

The model (1) will be estimated using maximum likelihood (ML). Maximum likelihood estimates have desirable asymptotic properties (consistency, efficiency, and asymptotic normality) and are robust to distributional assumptions for realistic departures from the Poisson assumption (note that the multinomial distribution leads to identical estimates with the Poisson distribution). Furthermore, they are essentially unbiased even for a very small sample of flows (Sen and Smith 1995).

Under some mild conditions (Sen and Smith 1995), the ML estimate of ζ = (A (1),...,A (_i ), B (1),...,B (J ), θ1 ,...,θ Q )´ exists and that of θ is unique. The estimates of A and B are not unique; however, if one A (i ) or B (j ) is chosen to be an arbitrary positive number, the remaining A (i )'s and B (j )'s are unique under the previous mild conditions. It can be proved that the ML estimate of ζ results in the solution to the following system of equations

(3)

(4)

(5)

where the operator * indicates summation with respect to the subscript it replaces (e.g., T_i* = Σ_jT_ij, T_*j = Σ_iT_ij ). A number of standard numerical methods or more specialized procedures can be used to solve equations (3) through (5). The three procedures adapted for the paper are the Deming-Stefan-Furness (DSF) procedure, the linearized DSF procedure, and the Modified Scoring procedure. An account of the development of the three procedures is given in Sen and Smith (1995). For completeness of the presentation, some of the details for each procedure follow.

The DSF Procedure for Parameters A_i and B_j

The DSF procedure gives values of T_ij for any choice of (θ' ), and N_i* ,N_*j , with Σ_iN _i* = Σ_jN_*j . Generally speaking, this procedure adjusts the rows (columns) of a two-dimensional table in each even (odd) iteration. After choosing an initial value for the column balancing coefficient, , say, the DSF procedure iterates as follows (the index r denotes the iteration number):

(6)

(7)

where O_i = T_i* , Dj = T_*j , and F_ij is a function of the separation measures . Upon convergence (see Sen and Smith (1995) for a proof of convergence), the values of T_ij are given by

.

In this paper, we chose a fairly stringent criterion for convergence as follows

(8)

where δ = 10^–12 . The algorithm attained this criterion in less than 100 iterations.

The DSF procedure essentially expresses T_ij as a function of θ. These values of T_ij could then be used in equation (5) to solve for an updated value of θ. In general, however, as θ changes, equations (3) and (4) could be violated, unless these changes are very small. This is achieved by using a linearized version of the DSF procedure, called the LDSF procedure (Weber and Sen 1985), which is computationally very attractive.

The LDSF Procedure for T_ij

Let us assume that we have run the DSF procedure and obtained a good set of T_ij that solves equations (3) and (4) for any given θ. This means that O_i = T_i* = N_i* and D_j = T_*j = N_*j . Define by ΔΟ = (ΔΟ₁,...,ΔΟ_I)' and ΔD = (ΔD₁,...,ΔD_J)' to be small changes in the values of Ο = (Ο₁,...,Ο_J)' and D = (D₁,...,D_J)' , respectively. Also, let ΔF_ij , be a small change in F_ij = exp[ θ'C_ij] . It can be proved (Sen and Smith 1995) that the corresponding small change, ΔT_ij , in each T_ij , so that ΔT_i* = ΔΟ_i , and ΔT_*j = ΔD_j can be obtained by the LDSF procedure, which iterates as follows

(9)

(10)

for , and r = 0, 1, 2,..., with initial T_ij values given by (11)

A proof for the convergence of the procedure is given in Weber and Sen (1985).

Changes in T_ij as a Function of a Change in θ

For a small change Δθ in θ and a small change 0 so that ΔΟ = ΔD = 0, it can be proved (Yun and Sen 1994) that an approximation for the corresponding small change ΔT_ij for each T_ij for all and is given by equation 12.

 (12)

The 's are constants, and O_i = T_i* , Dj = T_*j . Therefore, if T_ij 's are known, the only unknown in equation (12) is the Δθ . The solution for the Δθ will be the topic of the next section.

Estimation of θ Using the Modified Scoring (MS) Procedure

So far, for an initial value for θ, we obtained T_ij (θ) by using the DSF procedure to solve equations (3) and (4). We then changed θ to θ + Δθ and computed, using the LDSF procedure, with Δ0 = ΔD = 0, the corresponding change, ΔTij(θ, Δθ) in T_ij (θ) as a function of Δθ.

We are ready now to insert the [T_ij(θ) + ΔT_ij(θ, Δθ)]'s into the left side of equation (5) and solve the resultant equation for Δθ. Inserting T_ij + ΔT_ij in place of T_ij in (5) and using (12), equations (3) and (4) would remain approximately satisfied, while obtaining the following system of Q linear equations with Q unknowns (the Δθ_q's):

(13)

This system of equations can be solved by any standard solution method such as Gaussian elimination.

The current solution for θ at iteration r is updated next using the formula

θ^r = θ^r-1 + Δθ^r-1 (14)

If the corrections Δθ^r-1 have become negligible, the values of θ have been stabilized and the MS procedure terminates. Otherwise, new T_ij 's are obtained from the DSF procedure and the MS procedure continues. There is no guarantee that the MS procedure always converges (Sen and Smith 1995), although our computational experience is positive.

Goodness of Fit

Under the previous assumption that observations N_ij are independently Poisson distributed, the (Pearson) X² statistic,

(15)

where is an estimate of T_ij , is an appropriate measure of the overall fit of a model. Moreover, when is obtained using maximum likelihood, equation (15) has a X² distribution with df = ij – i –J – K + 1 degrees of freedom (Bishop et al. 1975; Rao 1973).

If , then X² = Z² , where

(16)

Since E (N_ij ) = T_ij and because N_ij have the Poisson distribution,

var(N_ij) = E(N_ij – T_ij)² = T_ij (17)

Therefore, E (Z²) = IJ , where I is the number of origin zones i, and J the number of destinations j . Equivalently, E (Z² /IJ ) = 1. Thus, the so-called "X² -ratio," X² /df , has an expectation that is asymptotically 1. It can be shown (Sen and Smith 1995) that the variance of the X² -ratio is

(18)

Hence, if T_ij 's are bounded away from zero (which is the case in exponential gravity models with finite parameters θ), the variance of Z²/IJ→0, as IJ→∞. It follows that when and T_ij 's are bounded away from zero, the variance of X²/df→0.

In practical applications, since the Poisson assumption seldom holds perfectly (as is the case here, where every pound or dollar value of shipment does not travel independently of each other pound or dollar value), an X² ratio less than 2 is a good indication that the gravity model fits the data well (Sen and Smith 1995).

Covariance of Maximum Likelihood Estimates

Covariance of `s

Let small case letters stand for the logarithms of corresponding capital letters (e.g., t_ij = log [T_ij ], a (_i ) = log [A (_i )], b (j ) = log [B (j )]. The model (1) may be written as

(19)

Let M denote the coefficient matrix of the right side of the system of equations (19). The matrix M is not of full rank (Sen and Smith 1995). However, the matrix M₍₂₎ obtained by deleting one of the first I + J columns of M is of full rank and has dimension IJ × (I + J + Q – 1). Let diag( · ) stand for a diagonal matrix, the diagonal elements of which are given within the parentheses. Then compute the matrix M'₍₂₎ · diag(T) · M₍₂₎ from the equation

(20)

where

(21)

(22)

U3 = ((upq )) with , with , with , V₁ = diag(T_1*,...,T_I*) , V² = ((T_ij)) and V₃ = diag(T_*1,...,T_{*j – 1}). Notice that the subscript j in each of the matrices above goes only up to J – 1.

Matrix M'₍₂₎ · diag(T) · M₍₂₎ , a square matrix of dimension (I + J + Q – 1), is the covariance matrix of M'₍₂₎ · diag(N). This is because the N_ij's have independent Poisson distributions and the covariance matrix Cov(N ) of N is diag(T ). It can be shown (Sen and Smith 1995) that the covariance matrix of is

Φ^-1 · M'₍₂₎ · diag(T) · M₍₂₎ · (Φ^-1)' (23)

where

Φ = M'₍₂₎ · diag(T) · M₍₂₎ · diag(1/A(1),...,1/A(I), 1/B(1),...,1/B(J – 1), 1,...,1) (24)

and

Φ^-1 = diag(1/A(1),...,1/A(I), 1/B(1),..., 1/B(J – 1), 1,...,1) · (M'₍₂₎ · diag(T) · M₍₂₎)^-1 (25)

Notice that using equations (23) through (25), another expression for the covariance matrix of can be written by as follows

diag(1/A(1),...,1/A(I), 1/B(1),..., 1/B(J–1),1,...1) ·  

(M'₍₂₎ · diagT · M₍₂₎)^-1 · diag(1/A(1),...,1/A(I), 1/B(1),..., 1/B(J–1), 1,...,1) (26)

The bottom right Q × Q submatrix of matrix (26) is the estimated covariance matrix of . The bottom right Q × Q submatrix of the inverse of (20) is (Rao 1973)

(27)

From equation (26), it follows that equation (27) is the covariance matrix of .

Covariance of _ij 's

Having obtained Cov(A ,B ,θ) from equation (23) or using equation (26), and since B (J ) is set equal to a constant, and its variance and covariances involving it are zeros, it can readily be seen that the covariance matrix of , denoted by the symbol Cov() , is

Cov() = Ψ₁ · Cov(A,B,θ) · Ψ'₁ (28)

with

Ψ₁ = diagT · M₍₂₎(1/A(1),...,1/A(I), 1/B(1),..., 1/B(J), 1,...,1) (29)

and M the coefficient matrix of the right side of the system of equations (19).

Short-Term Forecasting

The shipment of goods is affected by a multitude of factors (USDOT 1996, table 2.1). The weight and value of commodities shipped, characteristics of immediate concern in this study in particular, may be affected by the economy as a whole, globalization of business, international trade agreements, just-in-time inventory practices (weight only), packaging materials (weight only), economic regulation/deregulation, publicly provided infrastructure (weight only), user charges and other taxes, changes in truck size and weight limits (weight only), and technological advances. The discussion below is based on the assumption that, in the short term (e.g., three to five years) the compound effect of these factors on the size of the weight and value characteristics of commodity shipments is consistent with Poisson randomness.

Let the random variable be a future observation¹ of the flow from i to j . In the short run, under the assumption that the separation configuration will not change during the forecast period, we may conjecture that and N_ij will be highly (serially) correlated. Thus, we may then argue that the variance of the difference of future and present observations will be smaller than the variance of future observations alone. Hence,

var( − N_ij) = var() + var(N_ij) − 2Cov(N_ij) (30)

where var(·) stands for the "variance of" and Cov(·) for the "covariance of," will be small. In such cases, especially if N_ij 's are known, Sen and Smith (1995, chapter 5) suggest that it would be preferable to predict ( − N_ij)'s and add the predictions to the N_ij 's. A possible set of predictions for ( − N_ij)'s are ()'s. Thus, if the future is not too far off,

N_ij + () (31)

may yield a better prediction than . The computation of can be made easily using the LDSF procedure described earlier.

Other Forecasts

Four separate cases are discussed here. The first case assumes that would be available exogenously, θ^ƒ would be the estimate from the base period assumed to remain unchanged into the forecast period, and the estimates A^f (i ) and B^f (j ) could be base period estimates assumed to remain unchanged into the forecast period, or one or both of them could be exogenous. Then, T_ij for the forecast period may be estimated by

(32)

The calculation of Cov( ) requires the computation of the covariance matrices of the A^f (i )'s, B^f (j )'s, and θ^ƒ. For those estimates that remain unchanged from the base period, the appropriate covariance matrix is the one for the base period. For estimates obtained exogenously, the covariance matrix needs to be supplied exogenously. Covariances of estimates obtained from different independent data are usually assumed to be zeros.

The second case assumes that and θ^ƒ would be available as before and the estimates B^f (j ) and could be base period estimates assumed to remain unchanged into the forecast period or could be obtained exogenously. Then, T_ij for the forecast period may be estimated by

(33)

The second case assumes that and θ^ƒ would be available as before and the estimates A^f (i ) and could be base period estimates assumed to remain unchanged into the forecast period or could be obtained exogenously. Then, T_ij for the forecast period may be estimated by

(34)

The computation of the covariance matrices for the second and third cases is similar to the first case once the appropriate Jacobians with respect to A^f (i )'s, B^f (j )'s and θ^ƒ have been obtained. Recall that equation (29) is the respective Jacobian for the first case.

A fourth case, arising when , θ^ƒ, , and are available, requires the use of the DSF procedure to generate the A^ƒ (i ) and B^ƒ (j ). Then, T_ij for the forecast period may be estimated by equation (32). The computation of the covariance matrices for this case is treated in detail by Sen and Smith (1995, p. 440).

It is interesting to observe at this point that origin-destination travel times and/or costs (which depend on routes chosen and traffic congestion on the routes) usually become available only after the gravity model has been applied to forecast flows. This difficulty affects only the forecasting phase (not the estimation phase, when we use base period data for which travel times can be observed or computed using available information) and can be alleviated by combining the distribution and assignment stages of the process into a single stage.

EMPIRICAL ANALYSIS

Data Issues

Freight Shipments

Oak Ridge National Laboratory (ORNL) provided two sets of origin-destination flow data for freight shipment weight and value between origin seaports/border ports of entry i and destination states j from the following sources:

• The Transborder Surface Freight Database. This database (http://www.bts.gov/transborder/prod.html) is distributed by the Bureau of Transportation Statistics and contains freight flow data by commodity type and surface mode of transportation (rail, truck, pipeline, or mail) for U.S. exports to and imports from Canada and Mexico.
• The Port Import Export Reporting Service (PIERS) Database. This commercial database (http://www.piers.com/about/default.asp), offered by the Journal of Commerce , offers statistics on global cargo movements transiting seaports in the United States, Mexico, and South America to companies around the globe.

Two matrices, N_ij , were developed: one for shipment weight and one for shipment value for each of the two databases with dimensions 128 × 50 and 144 × 50 for the Transborder and PIERS databases, respectively.

The limitations of these two data sources are well documented (Meyburg and Mbwana 2002). A caveat for the Transborder Database is that it is a customs, not a transportation, dataset, resulting in inconsistencies between U.S. and Canadian data and issues related to the accuracy of transshipment data.² One caveat for using PIERS is that reported origins and destinations may be billing addresses rather than shipment points.

Separation Measures

In identifying specific types of spatial separation that tend to impede or enhance the likelihood of interactions between points of entry and destination sites, the most obvious type involves physical space, as exemplified by travel distance and travel time, which are quantifiable in terms of meaningful units of measurement. ORNL provided two sets with separate impedance matrices, , k = 1, 2 from origin point i to destination county c .

The first dataset included impedances between origin seaports and destination counties, the second set between origin border ports and destination counties. Both measures were computed from ORNL's National Highway Network (NHN). The first of the two measures is the route (network) distance in miles; the second computes a function of travel time for different functional classifications of highway segments and adds these time penalties together while tracing the previous network routes. The previous impedance matrices presented two problems: 1) more than 50% of the cells were empty because not all seaports/border ports of entry are connected to each county through the NHN; and 2) occasionally, freight shipments moved between origin-destination pairs with a missing connection (separation measure).

In order to deal with both problems, we devised and implemented the following stepwise procedure for the distance matrix (the first separation measure, ).

1. Centroid coordinates were computed for all counties nationwide.

2. Spherical distances in miles from each county to every other county.

3. For every destination county with a missing connection to an origin point (seaport/border port of entry), the closest county with an existing connection to an origin point was estimated based on the previous county-to-county distances.

4. The missing connection from an origin point to a destination county was finally computed to be the distance between the same point of origin and the closest destination county augmented by 130% the (airline) distance between the two destination counties (as a proxy to the actual road distance). This simplification was made under the assumption that a missing connection would imply that the destination county is off the NHN and thus would require additional time to be accessed from its closest county on the NHN.

The above procedure resulted in the construction of a synthesized separation measure comprising both route distance (the ORNL estimate) and multiples of airline distance (our estimate). It is important to take a closer look at those two components of the (re-estimated) first separation measure .

The airline distance between origin point i and destination county c , even augmented by 30%, is only an approximation of the actual miles traveled and represents a surrogate for the complex set of factors that express the difficulty of overcoming separation. On the other hand, the route distance for the same origin-destination pair is the total over-the-road distance on a realistic route. As a measure of separation, the route distance provides better accuracy for the eastern half of the limited access highway system of the United States, which has a higher level of complexity resulting in less circuitous routes than in the western half.

Whereas the distance measure developed above is a measure that depends more or less on the physical characteristics of the road link, the impedance measure, a function of link travel time as estimated by ORNL, depends on the special roadway type or certain conditions encountered on the link. ORNL's estimating procedure could not be replicated in this paper for the missing impedances. A surrogate value was computed based on the average speed obtained by ORNL's distance and time estimates, and the synthesized distances.

The procedure classified ORNL's distance and time estimates into 10 deciles. The stratum average speed was computed by the ratio of the average distance and time for that stratum (table 1). Each missing time impedance was then estimated by the ratio of its corresponding distance estimate and the average speed of the stratum to which the distance estimate belongs.

Clearly, this procedure imposes some inconsistency in the imputation of travel times, because the average speed estimate is based on decile groupings of the initial estimates. The consistency of the procedure could be improved if in a second iteration the decile groupings are based on all distances (initial and imputed), but I did not attempt this for this paper.

The previous steps resulted in two re-estimated sets of separation measures: 1) a distance matrix and a travel time matrix between origin seaports and destination counties; and 2) a distance matrix and a travel time matrix between origin points of entry and destination counties. The final step in preparing the separation measures for model estimation was to estimate these measures between origin points and destination states rather than counties. This was made possible simply by computing for each origin point i and all destination counties in state j the average distance and travel time. That is,

, ,

and

,

where n is the cardinality of Cj . Note that had we avoided the intermediate step of origin-to-destination county distance estimation, say, by computing distances to state geographical centers, origin-destination pairs with both ends in the same state would have been indistinguishable in terms of their separation.

Other Relevant Variables

As mentioned earlier, the weight and value of commodities shipped may be affected by a nexus of factors that are either origin-specific, destination-specific, or solely dependent on separation between an origin and a destination. Although study resources did not permit the collection of pertinent data on these factors, it could be interesting to observe how such factors can be accommodated in the proposed theoretical framework for future reference.

For origin-specific and destination-specific factors, Sen and Smith (1995) propose the following exploratory analysis: estimate A_i 's and Bj 's first using the maximum likelihood procedure described earlier, and then use these estimates as the dependent variable values in a model fitting procedure. Thus, the estimates of the A_i 's can be associated with the origin-specific factors, while the estimates of the Bj 's can be associated with the destination-specific factors.

Factors that can be solely attributed to the separation between an origin and destination can be accommodated in a deterrence function, an example of which is indeed the most general exponential function in equation (2). Openshaw and Connolly (1977) have compiled a list of possible functions and empirically compared several of them.

Results

The only decision necessary to apply the procedures described above (i.e., DSF, LDSF, and MS) is the choice of the flow unit. In the case of passenger transportation, a flow unit of 1 (consistent with a Poisson or multinomial distribution assumption) would be reasonable. In the case of freight shipments of goods, a basic unit of flow would appear to be a trainload (for shipments by rail) or a truckload (for shipments by truck). In the absence of mode-specific information as well as information related to the variation in modal size, we experimented with different values and determined an "optimal" (with regard to providing the best model fit) basic unit of flow of 100,000 pounds or dollars. This is not surprising given that the bulk of flows are long-distance shipments usually performed by large³ trucks or rail (on a limited scale for the particular data) with an average shipment value of $2.03 per pound for the Transborder data and $1.39 per pound for the PIERS data. Interestingly, our results appeared to be quite insensitive to the choice of the flow unit, primarily because the flows were inordinately large. This is in agreement with previous work (Sen and Pruthi 1983).

The procedures described above were run to a tight convergence. For each iteration of the MS procedure, the DSF procedure attained a 10E –12 convergence as defined by equation (8) in less than 100 iterations. The MS procedure itself attained a 10E –06 convergence as defined by the right-hand side of equation (13) in less than 20 iterations.

Parameter Estimates and Model Fit for PIERS Data

Separate gravity models were estimated for both seaport-to-state (seaport data) and border port-to-state (transborder) weight and value flows. In particular, for the seaport data, several model specifications were tested based on transformations of the distance and time separation measures, and , respectively, after a careful examination of residuals that showed the presence of outliers.

In the analysis of residuals, two difficulties were addressed. The first difficulty is the unequal variances of residuals, a consequence of the Poisson distribution. This problem was handled by considering, as Cochran suggests (Rao 1973, p. 393), instead of the plain residuals , the components

(35)

The second difficulty stems from the very large number of residuals (the number is I × J , the number of all origin-destination pairs). To tackle this concern, we followed a procedure described in Sen and Smith (1995, p. 456) consisting of drawing normal (rankit) plots of residuals. If the residuals are normal, such a plot would lie approximately on a straight line. Points that deviate sharply near the ends were tagged as outliers. Unusual jumps in the plots or other strongly nonlinear shapes could signal the need for transformations of the 's or for additional 's. There are clearly other ways of examining residuals. It is important to note here that the effort to conduct residual analysis in the gravity model case should not be seen as less than that for linear models where residual analysis is routinely carried out.

Indeed, the examination of residuals signaled the need for transformation of the 's. We found the square root transformation provided adequate fit for both weight flow and value flow estimates (table 2). The model provides an excellent fit for the data, because the Chi-square statistic hovers around its expected value of 1. In addition, the cell-to-cell weight flow estimates correlate very well with the data (Pearson correlation coefficient, r = 0.89, p < 0.01. The cell-to-cell value flow estimates also correlate very well with the data (r = 0.98, p < 0.01). Moreover, the weight and value length distributions in figure 1 (1% corresponds to 2,292,872,822 pounds) and figure 2 (1% corresponds to $3,198,123,011), respectively, seem to corroborate the previous results.

Of interest in table 2 is the appearance of with a positive value. This is due to collinearity between the used impedances, distance and travel time, a phenomenon that is known to adversely affect the sign of parameter estimates. Dropping one of the impedances would bias the parameter estimates left in the model. Given the robustness of the maximum likelihood procedure above in collinearity situations, all available impedance measures were retained (Sen and Smith 1995, chapter 5). After all, the sign of would have changed to a negative value had we reparameterized the model and considered instead of travel time as the first impedance measure, the difference between distance and travel time (in appropriate units).

In the end, the model specifications for the weight and value flows, respectively, are

(36)

Parameter Estimates and Model Fit for Transborder Data

In the case of Transborder data, we conducted a similar investigation of residuals (as above with the PIERS data) and observed the presence of a few remaining outliers, despite using transformed 's. These outliers were mainly states receiving an unusually large share of shipments compared with other states. The discovery was made after squaring the residuals in equation (35) and adding them together for all origin points of entry for each destination state. Destination states with a large share of shipments had a much larger sum of squares of residuals than other destinations. This was understandable, because these large destination states attract more shipments from longer distances and suggest the use of an additional impedance variable

for the weight flows, and

for the value flows, where δ is 1 for destinations with more than a 10% share in weight (4% in value) and 0 otherwise. In effect, we sought a different parameter estimate for the square root/log of distance for destination states with large shares.

The previous use of the indicator (dummy) variable δ is typical in spatial analysis for the treatment of residuals (Sen and Smith 1995) and introduces the effects of spatial structure on flow patterns (Gensler and Meade 1988; Fotheringham and O'Kelly 1989; Lo 1991), which is not the primary intent of this paper. More complex measures of relative location have been developed and tested empirically (Boots and Kanaroglou 1988; Lowe and Sen 1996). In the end, the model specifications for the weight and value flows, respectively, are

(37)

The previous steps removed all remaining outliers. Parameter estimates and goodness-of-fit statistics are shown in table 3. The previous observation regarding the sign of applies. The model fits the data very well, as the Chi-square ratios for both weight and value remain under 2. In addition, the cell-to-cell weight flow estimates correlate very well with the data (Pearson correlation coefficient, r = 0.91, p < 0.01). The cell-to-cell value flow estimates also correlate very well with the data (r = 0.94, p < 0.01). Moreover, the weight and value length distributions in figure 3 (1% corresponds to 767,519,815 pounds) and figure 4 (1% corresponds to $1,565,313,371), respectively, seem to corroborate the previous results.

Computation of Covariance of Estimates for PIERS Data

The previous point estimates of the parameters and the flows were used in the methodology described earlier to obtain the covariance of these estimates. The procedure can accommodate any reasonable number of separation measures, , and the large number of origins and destinations typically encountered in practice.

The computational requirements of the procedure are no longer prohibitive. We were able to obtain the covariance matrix of a 144 × 50 flow matrix, a 7200 × 7200 matrix, in a little more than an hour on a Pentium III, 800 MHz, 512MB RAM laptop computer running a FORTRAN 77 compiler. Replacing T_ij 's by their estimates in equation (27), the covariance matrix of for the weight flows from the PIERS dataset was found to be

(38)

The correlation between (the parameter estimate for the distance measure) and (the parameter estimate for the travel time measure) is readily apparent from equation (38). The negative sign of the covariances should not be disconcerting. It shows that, if for some small shift in the observations, were to increase, would decrease to "compensate."

Similarly, the covariance matrix of for the value flows from the PIERS dataset was found to be

(39)

The correlation between and is similarly apparent from equation (39).

Computation of Covariance of Estimates for Transborder Data

The covariance matrices of from the Transborder dataset for the weight and value flows, respectively, were found to be

(40)

and

(41)

These last two covariance matrices in equations (40) and (41) show relatively high correlation between (the parameter estimate for the distance measure) and (the parameter estimate for the travel time measure) and relatively low correlation between (the parameter estimate for the transformed travel time measure) and or . The interpretation of the negative signs is the same as above.

For illustration purposes, an example of the usefulness of the computed covariance matrices for 's in computing confidence intervals is shown in table 4. The table shows a few N_ij 's, 's, percentage difference between N_ij and , and 90% confidence intervals constructed by adding and subtracting 1.65 times the standard error available from the variances.

CONCLUSIONS

The most characteristic and indeed more restrictive feature of the Poisson distribution is the equality between its mean and variance. Many types of spatial interaction phenomena exhibit restricted (variance less than their means) or extra variation (variance greater than their means). Yet, the framework described above is very robust. Its robustness stems from the many types of asymptotic results that establish the Poisson distribution as the unique limiting form for a wide range of interaction processes when the population size increases and the overall influence of each individual interaction decreases.

More specifically, this research demonstrates that, to the extent the behavior of each shipping firm can be thought of as a very small interaction process, the resultant interaction process (i.e., the realized pattern of origin-destination commodity shipments) can be thought of as the superposition of all these processes and thus characterized by gravity models. The limited empirical analysis in this paper confirmed this conjecture and thus establishes an extremely rich framework for future experimentation.

Future applications of the proposed framework could conduct additional exploratory analyses to determine which factors affecting the demand for freight shipments are origin-based, destination-based, or separation-based. Another important issue needing more attention in the future is the requirement that, when the proposed modeling framework is used for forecasting, the travel times and costs used as impedance measures should be consistent with those obtained during the traffic assignment.

In summary, this paper has provided evidence that the movement of freight shipments can now be estimated within a desired confidence level as a result of maximum likelihood estimation of Poisson gravity models. The freight transportation modeler has two procedures available for computing reliable information within a predetermined accuracy: one for computing freight flow estimates and one for computing covariance matrices. These procedures can accommodate any reasonable number of separation measures,, and the large number of origins and destinations typically encountered in practice.

ACKNOWLEDGMENTS

This paper reports on research conducted by the author for the Bureau of Transportation Statistics, U.S. Department of Transportation.

REFERENCES

Bishop, Y.M.M., S.E. Fienberg, and P.W. Holland. 1975. Discreet Multivariate Analysis: Theory and Practice. Cambridge, MA: MIT Press.

Boots, B. and P. Kanaroglou. 1988. Incorporating the Effects of Spatial Structure in Discrete Choice Models of Migration. Journal of Regional Science 28:495–507.

Cambridge Systematics. 1997. A Guidebook for Forecasting Freight Transportation Demand, NCHRP Report 388. Transportation Research Board, National Research Council.

Fotheringham, A. and M. O'Kelly. 1989. Spatial Interaction Models: Formulations and Application. Studies in Operational Regional Science. Boston, MA: Kluwer Academic Publishers.

Gensler, W. and M. Meade. 1988. Locational and Population Factors in Health Care-Seeking Behavior in Savannah, Georgia. Health Services Research 22:443–462.

Lo, L. 1991. Spatial Structure and Spatial Interaction: A Simulation Approach. Environment and Planning A 23:1279–1300.

Lowe, J. and A. Sen. 1996. Gravity Model Applications in Health Planning Analysis of an Urban Hospital Market. Journal of Regional Science 36:(3)437–462.

Memmot, F.W. 1983. Application of Statewide Freight Demand Forecasting Techniques, NCHRP Report 260. Transportation Research Board, National Research Council.

Metaxatos, P. 2002. International Trade Traffic Study (Section 5115): Measuring the Ton-Miles and Value-Miles of International Trade Traffic Carried by Highway for Each State: Estimation of Origin-Destination Flows and Confidence Intervals, Final Report, prepared for the U.S. Department of Transportation, Bureau of Transportation Statistics. March.

Meyburg, A.H. and J.R. Mbwana, eds. 2002. Conference Synthesis: Data Needs in the Changing World of Logistics and Freight Transportation. New York Department of Transportation, Albany, NY. Available at http://www.dot.state.ny.us/ttss/conference/synthesis.pdf.

Openshaw, S. and C.J. Connolly. 1977. Empirical Derived Deterrence Functions for Maximum Performance Spatial Interaction Models. Environment and Planning A 9:1067–1079.

Rao, C.R. 1973. Linear Statistical Inference and Its Application. New York, NY: Wiley.

Sen, A. and T.E. Smith. 1995. Gravity Models of Spatial Interaction Behavior. New York, NY: Springer-Verlag

Sen, A. and R.K. Pruthi. 1983. Least Squares Calibration of the Gravity Model When Intrazonal Flows are Unknown. Environment and Planning A 15:1545–1550.

U.S. Department of Transportation (USDOT). 1996. Quick Response Freight Manual, DOT-T-97-10. Washington, DC. September 1996.

Weber, J.S. and A. Sen. 1985. On the Sensitivity of Gravity Model Forecasts. Journal of Regional Science 25:317–336.

Yun S. and A. Sen. 1994. Computation of Maximum Likelihood Estimates of Gravity Model Parameters. Journal of Regional Science 34:199–216.

Author address: Paul Metaxatos, University of Illinois at Chicago, 412 South Peoria Street, Suite 340, Chicago, IL 60607. Email: pavlos@uic.edu.

KEYWORDS: freight origin-destination flow estimation, covariance of estimates, gravity model.

1. The superscript f denotes that the quantity is an estimate for the forecast period.

2. See http://www.bts.gov/ntda/tbscd/desc.html, as of January 2004, for more information.

3. Currently, the gross weight limit for a 6-axle combination truck is 80,000 pounds.