Methods And Algorithms
ESTIMATION METHODS
INTRODUCTION
The estimation methods owe largely to Fair (1984). These first two sections are drawn from this book that we highly recommend. Econometric models are generally nonlinear, simultaneous, and large. They also tend to have error terms that are serially correlated. The estimation methods of IODE try to handle these characteristics. The notations that will be used in this chapter are as follows. Write the model as:
where
\(y_t\) is an n-dimensional vector of endogenous variables,
\(x_t\) is a vector of predetermined variables,
\(a_i\) is a vector of unknown coefficients,
\(u_{it}\) is an error term.
Assume that the first \(m\) equations are stochastic, with the remaining \(u_{it} (i = m+1 \dots n)\) identically zero for all \(t\). Let \(J_t\) be the \(n \times n\) Jacobian matrix whose \(ij\) element is \(f_i/\partial y_{jt} (i, j = 1 \dots n)\). Also, let \(u_i\) be the T-dimensional vector \((u_{i1} \dots u_{iT} , ..., u_{m1} \dots u_{mT})^\prime\). Let \(a\) denote the k-dimensional vector \((a_1\prime \dots a_m\prime)\) of all the unknown coefficients. Finally, let \(G_t^\prime\) be the \(k_i \times T\) matrix whose column \(t\) is \(\partial f_i(y_t, x_t, a_i)/\partial a_i\), where \(k_i\) is the dimension of \(a_i\), and let \(G^\prime\) be the \(k \times m.T\) matrix,
where \(k = \sum^m_{i=1} k_i\). These vectors and matrices will be used in the following section.
SERIAL CORRELATION AND DYNAMIC ADJUSTMENT MODEL
A convenient way of dealing with serially correlated error terms is to treat the serial correlation coefficients as structural coefficients and to transform the equations into equations with serially uncorrelated error terms. This introduces nonlinear restrictions on the coefficients, but otherwise the equations are like any others with serially uncorrelated coefficients. Consider the ith equation of (eq_1) and assume that uit, is first-order serially correlated:
where \(e_{it}\) is not serially correlated. Lagging (eq_1) one period, multiplying by \(\rho_i\), and subtracting the resulting expression from (eq_1) yields:
The equation includes now \(y_{t-1}\) and \(x_{t-1}\) as variables, and \(\rho_i\) as a coefficient, but it is no more general than (eq_1), and thus one can deal directly with (eq_4) under the assumption that serial correlation has been eliminated through transformation. The same procedure can handle serial correlation of higher orders. With respect to testing for serial correlation, it is well known that the Durbin-Watson (DW) test is biased toward accepting the null hypothesis of no serial correlation if there is a lagged dependent variable in the equation.
ESTIMATION TECHNIQUES
Ordinary Least squares (OLS)
Two-Stage Least Squares (2SLS)
Three-Stage Least Squores (3SLS)
Joint estimation of a non simultaneous system
ORDINARY LEAST SQUARES (OLS)
The OLS technique is a special case of the 2SLS technique, where Di in (eq_5) and (eq_6) below is the identity matrix. It is thus unnecessary to consider this technique separately from the 2SLS technique.
TWO-STAGE LEAST SQUARES (2SLS)
2SLS estimates of \(a_i\) (say \(\hat{a}_i\)) are obtained by minimizing
with respect to \(a_i\), where \(Z_i\) is a \(T \times K_i\) matrix of predetermined variables. \(Z_i\) and \(K_i\) can differ from equation to equation. An estimation of the covariance matrix of \(\hat{a}_i\) (say \(\hat{V}_{2ii}\)) is:
where \(\hat{G}_i\) is \(G_i\) evaluated at \(\hat{a}_i\), \(\hat{\sigma}_{ii} = T_i^{-1} \sum_{t=1}^T u_{it} u_{it}\), and \(u_{it} = f_i(y_{t-1}, x_{t-1}, \hat{a}_i)\).
The 2SLS estimator in this form is presented in Amemiya(L974). It handles the case of nonlinearity in both variables and coefficients. In earlier work, Kelejian (1971) considered the case of nonlinearity in variables only. Bierens (1981,p. 106) has pointed out that Amemiya s proof of consistency of this estimator is valid only in the case of linearity in coefficients, that is, only in Kelejian’s case. Bierens supplies a proof of consistency and asymptotic normality in the general case.
THREE-STAGE LEAST SQUORES (3SLS)
3SLS of \(a\) (say \(\hat{a}\)) are obtained by minimizing
with respect to \(a\), where \(\hat{S}\) is a consistent estimate of \(S\) and \(Z\) is a \(T \times K\) matrix of predetermined variables. As estimate of the covarialce matrix of \(\hat{a}\) (say \(\hat{V}_3\)) is
where \(\hat{G}\) is \(G\) evaluated at \(\hat{a}\). \(S\) is usually estimated from the 2SLS estimated residuals. This estimator is presented in Jorgenson and Laffont (1974), and it is further discussed in Amemiya(1977). Both prove consistency and asymptotic normality of 3SLS.
JOINT ESTIMATION OF A NON SIMULTANEOUS SYSTEM
Two aims can lead to jointly estimate several equations : equality restrictions on coefficients of different equations, and the presumed existence of contemporaneous correlation of the residuals (Seemingly unrelated equations: SURE). In this case, the joint estimation of the equations increases the efficiency of the estimators. Algorithmically, the method used is the Zellner method which consists in minimizing a weighted sum of the residuals. The weighting matrix is made of the contemporaneous covariance matrix of the residuals estimated with the residuals of a first step OLS on the equations. The Zellner can then simply be considered as a special case of the 3SLS method where the matrix of instruments is an identity matrix. The matrix \(D\) above is identically equal to the matrix \(\hat{S}_1 \otimes I\).
FULL INFORMATION MAXIMUM LIKELIHOOD (FIML)
Under the assumption that \((u_{1t} \dots u_{mt})\) is independently and identically distributed as multivariate N(O,S), the density function for one observation is
where \(S^* = S^{-1}\) and \(s^*_{ij}\) is the \(ij\) element of \(S^*\). The Jacobian \(J_t\) has been defined before. The likelihood function of t he sample \(t = 1 \dots T\) is
and the log of \(L^*\) is
Since log \(L^*\) is a monotonic function of \(L^*\), maximizing \(log ~ L^*\) is equivalent to maximizing \(L^*\). The problem of maximizing \(log ~ L^*\) can be broken into two parts: the first is to maximize \(log L^*\) with respect to the elements of \(S^*\), and the second is to substitute the resulting expression for \(S^*\) into (eq_10) and to maximize this concentrated likelihood function with respect to \(a\). The derivative of \(log ~ L^*\) with respect to \(s^*_{ij}\) is
where \(s_{ij}\) is the \(ij\) element of \((S^*)^{-1}\). This derivative uses the fact that \(? = a^*_{ij}\) for a matrix \(A\). Setting (eq_11) equal to zero and solving for \(s^*_{ij}\) yields
since \(?\) and therefore \(?\). Substituting (eq_12) into (eq_10) yields
The \(T\frac{m}{2}\) term comes from the fact that
The first and last terms on the RHS of (eq_13) are constants, and thus the expression to be maximized with respect to a consists of just the middle two terms.
Since \(log|S^*| = log|S^{-1}| = -log|S|\), the function to be maximized can be written
where, as noted earlier, the \(ij\) element of \(S\) is \(s_{ij}\). FIML estimates of a are thus obtained by maximazing \(L\) with respect to \(a\). An estimate of the covariance matrix of these estimates (say \(\hat{V}_4\)) is:
where the derivatives are evaluated at the optimum.
Phillips (1982) has pointed out that Amemiya’s pr oof of consistency and asymptotic efficiency (1977) is based on an incorrect lemma. This is corrected in a later paper (Amemiya Ig82). Amemiya’s article (1977), as corrected, shows that in the nonlinear case FIML is asymptotically more efficient than 3SI^S under the assumption of normality. In the linear case FIML is consistent even if the error terms are not normally distributed, where “FIML” means the full information maximum likelihood estimator derived under the assumption of normality. In the nonlinear case this is not in general true, although it sometimes is. Phillips (1982) presents an example of a non linear model for which FIML is consistent for a wide class of error distributions.
COMPUTATIONAL METHOD FOR LEAST SQUARES TECHNIQUES
LINEAR AND NONLINEAR LEAST SQUARES
The least squares estimation routine has been de signed to be very general and involve all the least squares methods: linear and nonlinear, single or system of equations, uses of instrumental variables or not, with or without correction for contemporaneous autocorrelation of residuals. The core of the computation is made of an iterative process which will be shown with the help of a little example.
Assume that the following equation is to be estimated:
where \(y\) is a \(T \times 1\) vector of observations of the dependent variable, \(X\) is the \(T \times Q\) matrix of observations of the independent variables. The \(Q \times 1\) vector a of coefficient are to be estimated. \(e\) is a \(T \times 1\) vector of residuals, normally distributed with mean \(0\) and standard deviation \(s\). It is also assumed that \(E(ee^\prime) = s^2 I\), \(I\) being the \(T \times T\) identity matrix. The function f is not necessarily linear in the variables and in the coefficients.The estimator of \(a\) is the vector \(\hat{a}\) which minimizes the sum of squares of the residuals of the equation (eq_16):
Let \(S=e^\prime e\) be the sum of squares of the residuals. Equation (eq_17) can be linearized by a Taylor expansion limited to the first term:
where \(v=u(\hat{a}_0) + e\), and \(u(\hat{a}_0)\) represents the terms neglected in the Taylor expansion, \(\hat{a}_0\) is a starting value for \(\hat{a}\), and \(y_0\) is defined by \(y_0 = f(X, \hat{a}_0)\). Equation (eq_18) is linear in \(\Delta \hat{a}_1 = (\hat{a} - \hat{a}_0)\). These parameters can be estimated by OLS, i.e. by
A new value for the coefficients can then be computed
This operation will be iterated until \(\Delta \hat{a}_k\), with \(k\) the iteration number, becomes sufficiently close to zero. At this stage, \(\hat{a}_k = \hat{a}\), and then \(y_k = f(X, \hat{a})\), \(u_k = 0\) and \(S = e^\prime e\).
The iterative process could well not converge, although it can be proved that, if at each iteration \(\Delta \hat{a}_k\) is weighted by \(\lambda\), with \(|\lambda| < 1\), so that \(\hat{a}_{k+1} = \hat{a}_k + \lambda \hat{a}_k\), the process must converge for sufficiently small values of \(\lambda\), provided that there is a minimum.
It is very easy to approximate the derivative of the f function, this is done numerically. In the linear case, the approximation is exact and gives the matrix of regressor observations. The iterative process reaches a solution within two iterations. In the nonlinear case, the matrix of derivatives is no more the matrix of observations, and a progressive correction process of the coefficients is repeated until the correction to apply is sufficiently small, i.e. smaller than a convergence threshold given by the user. At this point, the estimators minimize a weighted sum of the residuals of the equations, the \(S = S(\hat{a})\) function, the complexity of which is determined by the degree of non linearity of the system of equations. Nothing guaranties in very special non-linear cases that the process will eventually converge, or converge within the maximum number of iterations which was given by the user.
A weight (\(\lambda\) in the foregoing) can be introduced to accelerate convergence for each coefficient, these weights are called step relaxation parameters or smoothing parameters, by default their values are one, this value can be updated in the equation update screen. The step relaxation parameters are important for another reason: if the value of one or more relaxation parameters is set to zero, the initial value of the corresponding coefficients will be unchanged. This is a way to constraint a coefficient to a predetermined value.
The same iterative procedure is applied for a system of equations, \(y_k\) and \(u_k\) are no longer vectors of the \(T\) computed values of the endogenous variable and residuals, but are matrices of rank \(T\) times the number of equations of the system.
2SLS and 3SLS are incorporated in the procedure by first calculating a weighting matrix (a metric) using the specified instrumental variables, and then using this metric in the iterative computation of the least squares estimators of \(\Delta \hat{a}\). In this case, (eq_22) becomes
with \(D_1 = I \otimes Z (Z^\prime Z)^{-1} Z^\prime\).
3SLS and Zellner method require an additional iterative process: the first iterative process allows to compute an estimate of the covariance matrix of the contemporaneous residuals, the second iterative process re-estimates the coefficients by correcting the metric matrix. (eq_22) becomes:
with \(D_2 = S^{-1} \otimes Z (Z^\prime Z)^{-1} Z^\prime\) and \(S\) estimated during the first step. The least squares estimation program is very general. It can handle linear or nonlinear equations, single equations or a system of equations. The estimated coefficients have the BLUE properties if the model is linear, but only the asymptotic properties are ensured for non-linear models.
The estimated covariance matrix of the coefficients is exact for the linear models, but only asymptotically exact for non-linear models.
FIML
The algorithm used in IODE to calculate FIML estimates is a Gauss-Newton one. Numerical derivatives of the likelihood function are computed like in the least squares algorithm. The weighting matrix is the same as in the 3SLS procedure, implying that the user must specify a meaningful set of instrumental variables. (Moreover, the 3SLS estimates are strongly recommended as starting values.)
STANDARD STATISTICS
A limited number of statistics are given as standard output by the estimation program. Per coefficient: the value of the coefficient, the standard error of the coefficient, and the t-statistics. In the non-linear case the standard errors of the coefficients only represent a lower limit of the true standard errors which will be reached asymptotically. When the model is a linear single equation with a constant term, the total variance of the endogenous variable is equal to the variance explained by the regression plus the residual variance. The \(R2\) statistics represents the explained share of the total variance in percents. This statistics is corrected for the number of coefficients. Out of this limited case, the interpretation of this statistics is difficult and hazardous.
The standard error gives a more straightforward interpretation of the quality of the adjustment.
SIMULATION ALGORITHM
THE GAUSS-SEIDEL TECHNIQUE
Most macro-econometric models are solved using the Gauss-Seidel technique. It is a simple technique and in most cases works remarkably well. This technique is used in IODE to simulate models. The Gauss-Seidel technique is easy to describe by means of an example.
Assume that the model consists of three equations, and let \(x_{it}\) denote the vector of predetermined variables in equation \(i\). The model is as follows:
The technique requires that the equations be rewritten with each endogenous variable on the LHS of one equation. This is usually quite easy for macro-econometric models, since most equations have an obvious LHS variable.
If, say, the LHS variable for an estimated equation is \(ln(y_{2t}/y_{1t})\), then \(y_{2t}\) can be written on the LHS by taking exponents and multiplying the resulting expression by \(y_{1t}\).
The program can do that job for you since, for each equation, it requires the definition of the endogenous variable specified by the equation, provided that the endogenous variable appears only once in the equation.
The endogenous variable can appear either on the LHS or in the RHS. Simulation of a model requires moreover that the endogenous variable be specified only once in the full set of equations. When these transformations of the equations are done, the model can be rewritten as follows:
In order to solve the model, values of the coefficients and the error terms are needed. Given these values and given values of the predetermined variables, the solution proceeds as follows.
Initial values of the endogenous variables are guessed. These are either the actual values or extrapolations from the previous period according to some rules. IODE allow various strategies. Given these values the equations can be computed successively and produce a new set of values for the endogenous variables. Given this new set of values, the equations can again be solved to get another set. Convergence is reached if for each endogenous variables the values on successive iterations are within some prescribed tolerance level. There is no guarantee that the Gauss-Seidel method will converge.
The advantage of the technique, however, is that it can usually be made to converge (assuming an actual solution exists) with sufficient damping. By damping is meant the following. Let \(y_{k-1}\) denote the solution value of \(y\) for iteration \(k-1\) and let \(\hat{y}_k\) denote the value computed by solving the equation on iteration \(k\). Instead of using \(y_k\) as the solution value for iteration \(k\), one can adjust \(y_{k-1}\) only partway toward \(\hat{y}_k\). If \(\lambda\) is 1, there is no damping.
The question of what to use as a stopping rule is not as easy as it might sound. The stopping rule can either be in absolute or relative terms. In absolute terms it is
and, in relative, it is
where \(\epsilon_i\) is the tolerance criterion for variable \(i\). If damping is used, \(y_{it,k}\) must be replaced by \(\hat{y}\). A global criterion can also be used
The problem is to choose the value for \(\epsilon_1\). It is inconvenient to have to choose different values of the tolerance criterion for different variables, and one would like to use just one value of \(\epsilon\) throughout. This is not, however, a sensible procedure if the units of the variables differ and if the absolute criterion is used. The problem is lessened if the percentage criterion is used, but in this case one must be concerned with variables that can be close to zero or can be zero. The solution adopted in IODE is to compute a stopping rule which combines the three afore-mentioned criteria
EQUATION REORDERING
To increase the simulation speed, the model can be split into Connex Components. Each component will be reorganized by an heuristic algorithm to obtain a quasi sub-diagonal model.
Let E denote the incidence matrix:
\(E_{ij} = 1\) if the variable \(j\) belong to the equation, 0 otherwise
After decomposition, the incidence matrix \(E\) has the following form:
CONNEX COMPONENTS CONSTRUCTION
Let VARS be the list of endogenous variables
Take the first remaining endogenous variable in VARS : ENDO0
Compute SUCC = list of variables depending on ENDO0 (recursively)
Compute PRE = list of variables from which depends ENDO0 (recursively)
Compute \(CFC_i\) intersection between SUCC and PRE : new connex component
Let \(VARS = VARS \ CFC_i\)
If \(VARS \neq \varnothing\) go to 1, else stop
Reordering of the connex components
Construct the \(E^\prime\) matrix:
To make \(E^\prime\) lower triangular, apply the following algorithm (\(O\) is the new order of Connex Components):
\(O = \{\}\)
Choose \(i \neq 0\) so that \(E^\prime_{ij} = 0, \forall j \neq i\)
Set \(E^\prime_{ij} = 0, \forall j\)
Add \(i\) in \(O\)
Go to 2 if at least one \(CFC\) remains
Triangulation of a Connex Component
Each Connex Component (CFC) can be reorganized to minimize the computation time required to solve the corresponding system. The following algorithm reduces approximately the number of iterations by a factor 2. It is based on the fact that to maximize the information during an iteration, the extra-diagonal elements of \(E\) must be as close as possible to the diagonal.
compute the sum of extra-diagonal elements:
\(\forall i ~ \exists ~ CFC\) set equation \(i\) after equation \(k : max_{k>i}, k : E_{ki} = 1\)
Compute \(S1^\prime\), \(S2^\prime\) and \(S3^\prime\) as \(S1\), \(S2\) and \(S3\). If \(S1^\prime < S1\) and \(S2^\prime < S2\) and \(S3^\prime < S3\), set \(S1 = S1^\prime\), \(S2 = S2^\prime\), \(S3 = S3^\prime\) and go to 2 else stop
Goal seeking
Given \(n\) pair \(ENDO_i-EXO_i\) variables, the solution of the system of equations must now be computed for the new set of variables consisting of the initial group in which the ENDO variables for the given list have been replaced by EXO variables. If we where to use a NEWTON algorithm or an as simulated method, we would just have to solve the system for the new set of variables. In the case of the GAUSS algorithm however one must declare a new output variable for each equation. As an example suppose we do one endo-exo exchange. (for more pai rs we just repeat the process). The goal is to find a pair of equations, the first of which contains EXO and the endogenous variable of the second being ENDO. We therefore construct a structured tree as below where an arrow indicates that a variable belongs to the definition equation of another variable.
