Overview   Marginal effects measure the expected instantaneous change in the dependent variable as a function of a change in a certain explanatory variable while keeping all the other covariates constant. The marginal effect measurement is required to interpret the effect of the regressors on the dependent variable. This example illustrates the calculation of marginal effects by using the QLIM procedure in binary choice models and censored models. Binary Choice Model   The first data are the binary choice (0–1) data set used by Spector and Mazzeo (1980) to study the effectiveness of a new method in teaching Economics. The data are read by the following DATA step statements. data greenedata; input gpa 1-4 tuce 6-7 psi grade; datalines; 2.66 20 0 0 2.89 22 0 0 ... more lines ... 3.1 21 1 0 2.39 19 1 1 ; run; The dependent variable y is modeled as follows: where   is the conditional mean function, x is the vector of explanatory variables, and   is the error term. The conditional mean function is given by where F denotes a cumulative distribution function and β denotes the parameters. Therefore, Marginal effect is a measure of the instantaneous effect that a change in a particular explanatory variable has on the predicted probability of y, when the other covariates are kept fixed. They are obtained by computing the derivative of the conditional mean function with respect to x given by where   is the density function that corresponds to the cumulative function  . The marginal effects are nonlinear functions of the parameter estimates and levels of the explanatory variables. Hence, they generally cannot be inferred directly from parameter estimates. Marginal effects for distributions such as probit and logit can be computed with PROC QLIM by using the MARGINAL option in the OUTPUT statement. The following MODEL statement fits the model equation to the endogenous variable GRADE and the covariates GPA, TUCE, and PSI. Further, you can specify the discrete nature of the endogenous variable by using the DISCRETE option. The D=PROBIT option in the MODEL statement enables you to specify the probit distribution. In the OUTPUT statement, use the OUT = option coupled with the MARGINAL option to obtain the marginal effects in the data set.   proc qlim data=greenedata; model grade = gpa tuce psi / discrete(d=probit); output out=outme marginal; run; quit;   The following output in Figure 2.1 displays the probit model fit statistics.     Figure 2.1 Parameter Estimates The QLIM Procedure   Parameter Estimates Parameter DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 -7.452320 2.542536 -2.93 0.0034 gpa 1 1.625810 0.693869 2.34 0.0191 tuce 1 0.051729 0.083891 0.62 0.5375 psi 1 1.426332 0.595036 2.40 0.0165   The MARGINAL option in PROC QLIM evaluates marginal effects for each observation. In the output data set, OUTME, 'Meff_P2_covariate' and 'Meff_P1_covariate' refer to the marginal effect of 'covariate' on the probability of GRADE=1 and on the probability of GRADE=0, respectively. Hence, 'Meff_P2_gpa' is the marginal effect of GPA on the probability of GRADE=1. To evaluate the "average" or "overall" marginal effect, two approaches are frequently used. One approach is to compute the marginal effect at the sample means of the data. The other approach is to compute marginal effect at each observation and then to calculate the sample average of individual marginal effects to obtain the overall marginal effect. For large sample sizes, both the approaches yield similar results. However for smaller samples, averaging the individual marginal effects is preferred (Greene 1997, p. 876). This example uses PROC QLIM to compute the overall marginal effect by these two approaches. PROC QLIM outputs the marginal effects computed at each observation in the data set. Hence, in order to obtain overall marginal effect, you can use PROC MEANS to obtain the sample average of individual marginal effects:     proc means data=outme n mean; var Meff_P2_gpa Meff_P2_tuce Meff_P2_psi; title 'Average of the Individual Marginal Effects'; run; quit;   Figure 2.2 shows the overall marginal effects of PSI, TUCE, and GPA for the probit model fit.   Figure 2.2 Average of Individual Marginal Effects   Average of the Individual Marginal Effects The MEANS Procedure   Variable Label N Mean Meff_P2_gpa Meff_P2_tuce Meff_P2_psi Marginal effect of gpa on the probability of grade=2 Marginal effect of tuce on the probability of grade=2 Marginal effect of psi on the probability of grade=2 32 32 32 0.3607863 0.0114793 0.3165199   Marginal effects at sample means of covariates can be obtained by adding an observation with a missing value for GRADE (endogenous variable) and the sample means of the covariates to the original dataset GREENEDATA. The following DATA step statements add the observation to the original dataset where means of GPA, TUCE, and PSI are 3.117, 21.938, and 0.4375, respectively. Because of the missing value of GRADE, model fitting does not depend on this observation. However, marginal effects and predicted probability values are obtained.   data me_mean; input gpa tuce psi grade; datalines; 3.117 21.938 .4375 . ; run; data me_mean; set greenedata me_mean; run;   Figure 2.3 shows the results of the following statements, which use the MARGINAL option in PROC QLIM on the new data set ME_MEAN:   proc qlim data=me_mean; model grade = gpa tuce psi / discrete(d=probit); output out=outme1 marginal; run; quit;     Figure 2.3 Marginal Effects at Covariate Means   Average of the Individual Marginal Effects Obs gpa tuce psi grade Meff_P1_gpa Meff_P2_gpa Meff_P1_tuce Meff_P2_tuce Meff_P1_psi Meff_P2_psi 33 3.117 21.938 0.4375 . -0.53325 0.53325 -0.016967 0.016967 -0.46783 0.46783   There might be cases where some regressors are dummy variables. Since the MARGINAL option in PROC QLIM uses derivative calculation, marginal effect results might not be meaningful. In some cases, it might only help to provide a reasonable approximation. In such cases, analyzing the marginal effects due to dummy variables by taking the difference of estimated probabilities between the different levels of dummy covariates (PSI=1 and PSI=0) is more appropriate. Let x k  denote the dummy explanatory variable and x* denote the other covariates at their means. The effect due to a change of x k  on the predicted probabilities of y is     The PROBALL option in the OUTPUT statement enables you to obtain the predicted probabilities of discrete endogenous variables for all responses. In order to obtain the previously described difference (when the other covariates are fixed at their means), you can append two observations with 1 and 0 as the discrete regressor’s (PSI) value, means of the other covariates, and missing values for the response variable as shown. data me_diff; input gpa tuce psi grade; datalines; 3.117 21.938 1 . 3.117 21.938 0 . ; run; data me_diff; set greenedata me_diff; run; proc qlim data=me_diff; model grade = gpa tuce psi / discrete(d=probit); output out=outme2 marginal proball; run; quit; The marginal effect of PSI is computed as follows: The table obtained by using the PROBALL option in PROC QLIM produces two columns labeled 'Prob2_grade' and 'Prob1_grade' which are   and   respectively. You can obtain the predicted probabilities required for the above computation of marginal effects from the last two observations in column 'Prob2_grade'. As shown in Figure 2. 4, you can compute the marginal effect of PSI to be 0.46437(0.57005 - 0.10568). Hence, the overall effect of a unit change in PSI (by introducing the new teaching method with an exposure to PSI), on the probability of GRADE=1 (an increase in the student’s grade) is 0.46.     Figure 2.4 Estimated Probability by Using the PROBALL Option   Average of the Individual Marginal Effects Obs gpa tuce psi grade Meff_P1_gpa Meff_P2_gpa Meff_P1_tuce Meff_P2_tuce Meff_P1_psi Meff_P2_psi Prob1_grade Prob2_grade 33 3.117 21.938 1 . -0.63858 0.63858 -0.020318 0.020318 -0.56023 0.56023 0.42995 0.57005 34 3.117 21.938 0 . -0.29701 0.29701 -0.009450 0.009450 -0.26057 0.26057 0.89432 0.10568   Note that although for discrete covariates the difference calculation is usually employed to obtain marginal effects, in this example you can observe that both the derivative (Meff_P2_psi= 0.46783, from Figure 2.3) and difference calculations (0.46437, by taking the difference of probabilities under 'Prob2_grade' from Figure 2.4) produce similar results. However, this might not be true for other data sets. The marginal effects plot with respect to PSI on the Pr(GRADE = 1) is shown in Figure 2.5  using results from the probit model fit. The marginal effects of PSI on Pr(GRADE = 1)  are obtained as a function of the GPA, at the mean of TUCE. This allows better interpretation of results. As shown in Figure 2.5, you can conclude that the effect of PSI is much greater on those students who have a higher GPA.     Figure 2.5 Marginal Effect of PSI on Predicted Probabilities     Similar analysis with ordered multinomial, censored, and truncated data are supported by the QLIM procedure. An example with censored data is illustrated below.   Censored Regression Model   In this section an illustration of the marginal effect calculation for censored regression model is given. The censored regression model is described as follows:   where   is the latent variable and yi is the observed variable. The marginal effect in this case is as follows:   The Mroz (1987) data set taken from the 1976 panel study of income dynamics based on data obtained for the previous year 1975, is used for analysis. The dataset is used in the study of female labor supply. A tobit model is fitted to the data. The following variables are used in the model: WHRS = wife’s hours of work in 1975 LFP = dummy variable; LFP = 1 if woman worked in 1975; else = 0 KL6 = number of children less than 6 years old in household K618 = number of children between ages 6 and 18 in household WA = wife’s age WE = wife’s educational attainment, in years HHRS = husband’s hours worked in 1975 HA = husband’s age HE = husband’s educational attainment, in years HW = husband’s wage, in 1975 dollars FAMINC = family income, in 1975 dollars The following DATA step statements are used to read the data.   data mroz; infile datalines truncover; input LFP WHRS KL6 K618 WA WE WW RPWG HHRS HA HE HW FAMINC MTR WMED WFED UN CIT AX; datalines; 1 1610 1 0 32 12 3.3540 2.65 2708 34 12 4.0288 16310 .7215 12 7 5.0 0 14 1 1656 0 2 30 12 1.3889 2.65 2310 30 9 8.4416 21800 .6615 7 7 11.0 1 5 1 1980 1 3 35 12 4.5455 4.04 3072 40 12 3.5807 21040 .6915 12 7 5.0 0 15 ... more lines ... In this example, a tobit model is estimated. You can fit a logit or a probit model to the dummy response variable LFP, similar to the analysis shown earlier. You can append the MROZ data set with the mean values of the explanatory variables as shown in the following statements. /*last row has the means*/ data add; input LFP WHRS KL6 K618 WA WE HHRS HA HE HW FAMINC ; datalines; 1 740.5763612 0.237715803 1.353253652 42.53784861 12.28685259 2267.270916 45.12084993 12.49136786 7.482178752 23080.59495 ; run; data mrozexsub; set mroz(keep=LFP WHRS KL6 K618 WA WE HHRS HA HE HW FAMINC) add; run;   Figure 2.6 shows the marginal effect results for a tobit model fit.   proc qlim data=mrozexsub; model WHRS = KL6 K618 WA WE HHRS HA HE HW FAMINC; endogenous WHRS ~ censored(lb=0); output out=outtobit residual marginal; run; proc print data=outtobit(firstobs=754); var Meff_KL6 Meff_K618 Meff_WA Meff_WE Meff_HHRS Meff_HA Meff_HE Meff_HW Meff_FAMINC; run;   Figure 2.6 Marginal Effects for a Tobit Model Fit     Obs Meff_KL6 Meff_K618 Meff_WA Meff_WE Meff_HHRS Meff_HA Meff_HE Meff_HW Meff_FAMINC 754 -579.163 -63.9936 -17.2332 52.7183 -0.43955 -10.7785 -3.79460 -127.807 0.039525     Similar analysis using the probit and logit models yield results shown in Figure 2.7 and Figure 2.8 respectively.     Figure 2.7 Marginal Effects for a Probit Model Fit   Obs Meff_P1_KL6 Meff_P1_K618 Meff_P1_WA Meff_P1_WE Meff_P1_HHRS Meff_P1_HA Meff_P1_HE Meff_P1_HW Meff_P1_FAMINC 754 0.34706 0.020687 0.011396 -0.056156 .000255579 .007300417 .006522166 0.068965 -.000021673       Figure 2.8 Marginal Effects for a Logit Model Fit Obs Meff_P1_KL6 Meff_P1_K618 Meff_P1_WA Meff_P1_WE Meff_P1_HHRS Meff_P1_HA Meff_P1_HE Meff_P1_HW Meff_P1_FAMINC 754 0.35548 0.019558 0.011089 -0.058672 .000274938 .008184995 .006452848 0.074250 -.000023536     References   Greene, W. H. (1997), Econometric Analysis, Third edition, Prentice Hall, 339–350. Mroz, T.A. (1987), "The Sensitivity of an Empirical Model of Married Women’s Hours of Work to Economic and Statistical Assumptions," Econometrica, 55, 765–799. Spector, L. and Mazzeo, M. (1980), "Probit Analysis and Economic Education," Journal of Economic Education, 11, 37–44. SAS Institute Inc. (2008), SAS/ETS User’s Guide, Version 9.2, Cary, NC: SAS Institute Inc. ... View more

Overview   Econometricians are often interested in the finite sample efficiency of different estimators when the underlying assumptions of least squares break down. This specific example compares the efficiency of ordinary least squares (OLS) with that of the Yule-Walker method (often called Prais-Winsten) when disturbance terms are assumed to be first-order autocorrelated. When first-order autocorrelation exists in OLS residuals, the OLS estimator is unbiased but not efficient. The sampling variances are underestimated, causing inferences from t- and F-tests to be invalid. One way to compensate for the autocorrelated residuals is the Yule-Walker method.The Yule-Walker method estimates the autoregressive form of the error term and then estimates the coefficients via generalized least squares (GLS). One particular advantage of the Yule-Walker method is that it retains the information in the first observation. This is crucial in small-sample estimation. The gain in efficiency can be investigated by performing a small Monte Carlo experiment. The data are generated by the following AR(1) model:     where   and u t  are disturbances terms. Based on this model, this example compares the finite sample efficiency of OLS and the Yule-Walker method by drawing 1000 samples with 20 observations each. Analysis In order to compare the finite sample efficiency of OLS and the Yule-Walker method by simulation, this example generates 1000 samples of 20 observations each for the following model: This generation of samples is accomplished by using the DATA step to turn the model into simulated data. For this example, the value of  is set equal to 0.9. You can choose to repeat the experiment several times with other values of   to investigate the relative efficiency of the estimators in different error persistence environments. data sample; do samp = 1 to 1000; rho = 0.9; /* first error term */ eps = rho * rannor( 47392 ) + rannor( 82745 ); do x = 1 to 20; y = 2 + 5 * x + eps; eps=eps; eps = rho * eps + rannor( 32815 ); output; end; end;   OLS   Each sample is then estimated by OLS using the AUTOREG procedure to specify a linear regression of Y on X and a CONSTANT. The OUTEST= option creates the SAS data set EST0 that contains the parameter estimates from each of the samples. The NOPRINT option suppresses printed output. Estimates are computed for each sample by using the BY statement to specify BY-group processing.   proc autoreg data=sample outest=est0 noprint; by samp; model y = x ; run;   The distribution of the parameter estimates can be calculated with the UNIVARIATE procedure. The following commands rename the two variables from the SAS data set EST0 and calculate their empirical distributions.   data estbar0; set est0; rename x=xbar0; rename intercep=inter0; run; proc univariate data=estbar0; var xbar0 ; histogram xbar0 / normal(color=blue) cframe=ligr cfill=green; title 'Distribution of OLS Estimate'; run;   Below are the summary statistics for the coefficient of X. The mean is very close to 5, the true value, and the variance is 0.022508.     Distribution of OLS Estimate The UNIVARIATE Procedure Variable: xbar0 (Parameter Estimate for x) Moments N 1000 Sum Weights 1000 Mean 5.00099049 Sum Observations 5000.99049 Std Deviation 0.15002553 Variance 0.02250766 Skewness 0.02338237 Kurtosis -0.2000697 Uncorrected SS 25032.391 Corrected SS 22.4851517 Coeff Variation 2.9999163 Std Error Mean 0.00474422 Basic Statistical Measures Location Variability Mean 5.000990 Std Deviation 0.15003 Median 5.000203 Variance 0.02251 Mode . Range 0.86598     Interquartile Range 0.20464 Tests for Location: Mu0=0 Test Statistic p Value Student's t t 1054.122 Pr > |t| <.0001 Sign M 500 Pr >= |M| <.0001 Signed Rank S 250250 Pr >= |S| <.0001 Quantiles (Definition 5) Quantile Estimate 100% Max 5.45220 99% 5.36620 95% 5.24473 90% 5.19162 75% Q3 5.10437 50% Median 5.00020 25% Q1 4.89972 10% 4.81080 5% 4.75461 1% 4.65464 0% Min 4.58623       Yule-Walker   The Yule-Walker estimator is calculated similarly. It is chosen as the method of estimation by specifying the NLAG= option in the model statement in the AUTOREG procedure. Again, the distribution of the parameter estimates is calculated by using the UNIVARIATE procedure.   proc autoreg data=sample outest=est1 noprint; by samp; model y = x / nlag=1; run; data estbar1; set est1; rename x=xbar1; rename intercep=inter1; run; proc univariate data=estbar1; var xbar1 ; histogram xbar1 / normal(color=blue) cframe=ligr cfill=green; title 'Distribution of Yule-Walker Estimate'; run;     Like the OLS estimator, the mean of the Yule-Walker estimator is very close to 5. Notice, however, that the variance of the Yule-Walker estimator is 0.019559, which is less than that of the OLS estimator. Also, an examination of the quantiles reveals less dispersion in the Yule-Walker distribution as a whole.   Distribution of Yule-Walker Estimate The UNIVARIATE Procedure Variable: xbar1 (Parameter Estimate for x) Moments N 1000 Sum Weights 1000 Mean 4.99968093 Sum Observations 4999.68093 Std Deviation 0.13985238 Variance 0.01955869 Skewness 0.02177098 Kurtosis -0.1709983 Uncorrected SS 25016.3485 Corrected SS 19.5391284 Coeff Variation 2.79722602 Std Error Mean 0.00442252 Basic Statistical Measures Location Variability Mean 4.999681 Std Deviation 0.13985 Median 4.995825 Variance 0.01956 Mode . Range 0.83045     Interquartile Range 0.19028 Tests for Location: Mu0=0 Test Statistic p Value Student's t t 1130.505 Pr > |t| <.0001 Sign M 500 Pr >= |M| <.0001 Signed Rank S 250250 Pr >= |S| <.0001 Quantiles (Definition 5) Quantile Estimate 100% Max 5.42195 99% 5.34807 95% 5.21860 90% 5.17968 75% Q3 5.09476 50% Median 4.99582 25% Q1 4.90448 10% 4.81394 5% 4.77470 1% 4.67258 0% Min 4.59151     Output   The output shows that, for this experiment, the Yule-Walker estimator is slightly more efficient than the OLS estimator. The results of the standard deviations of the distribution of the parameter estimates are s2yw = 0.019559 and s2ols = 0.022508 (where s2 denotes the variance). The fact that s2yw < s2ols implies greater efficiency for the Yule-Walker method. An F-test can be computed to see if the variances are significantly different from each other. The F statistic is the ratio of the variances for the two estimators adjusted for degrees of freedom, which is 999 for both the numerator and the denominator. The variance of the OLS estimator, s2ols, is 0.022508, while the variance of the Yule-Walker estimator, s2yw, is 0.019559. Thus, the appropriate F statistic is given by s2ols / s2yw = 0.022508 / 0.019559 = 1.1507746. The p-value of this statistic for a one-sided test can be computed with the following DATA step data _null_; p = 1-probf( 1.1507746 , 999, 999 ); put p=; run;   which returns a value of 1 - Pr F(1.1507746, 999, 999) = 0.013. Thus, you would reject the null hypothesis that the efficiency of the two estimators is the same for most conventional levels of significance.   References   Greene, W.H. (1993), Econometric Analysis, Second Edition, New York: Macmillan Publishing Company. Harvey, A.C. and McAvinchey, I.D. (1979), ``On the Relative Efficiency of Various Estimators of Regression Models with Moving Average Disturbances," Proceedings of the Econometric Society European Meetings, ed. E. Charatsis, Athens, 1979. Johnston, J. (1972), Econometric Methods, Second Edition, New York: McGraw-Hill. SAS Institute Inc. (1993), SAS/ETS User's Guide, Version 6, Second Edition, Cary, NC: SAS Institute Inc. ... View more

Overview   The Efficient Method of Moments (EMM) is a simulation-based method of estimation that seeks to attain the efficiency of Maximum Likelihood (ML) while maintaining the flexibility of the Generalized Method of Moments (GMM.) This is done using the scores of a pseudo (or auxiliary) model as moment conditions in the GMM step. EMM is particularly useful when the ML approach is not feasible or computationally intensive, such as in models with dynamic latent variables. In this example, the GMM feature of PROC MODEL is used for implementing the EMM method to estimate a simple stochastic volatility model, where a GARCH(1,1) model serves as the auxiliary model.   Details   EMM, introduced by Bansal et al. (1993 and 1995) and Gallant and Tauchen (1996), is a variant of the Simulated Method of Moments (SMM) of Duffie and Singleton (1993). The idea is to match the efficiency of the ML estimation with the flexibility of the GMM procedure. ML itself can be interpreted as a method of moment procedure, where the score vector, the vector of derivatives of the log-likelihood function with respect to the parameters, provides the exactly identifying moment conditions. In a first stage, EMM employs an auxiliary (or pseudo) model, which closely matches the true model, for the ML procedure. In a second stage, the scores of the auxilliary model, computed at the ML estimates, provide the moment conditions in the SMM procedure. EMM is useful for estimating models when the computation of the likelihood function is either infeasible or cumbersome. The typical application involves models with dynamic latent variables, where only one part of the state vector is observable and computation of the likelihood involves integrating out the unobserved part. The dimension of the integral is equivalent to the sample size, so that direct evaluation of the likelihood is not practical. An example of such a model, from Epidemiology, is the SEIR model that describes the parts of a population susceptible, exposed, infected, and recovered from a disease in terms of a system of differential equations, whereas usually the data report only those infected (Olsen and Schaffer 1990). Other examples include stochastic volatility models from Finance, where the instantaneous volatility is unobserved and only the security price can be measured (Gallant and Tauchen 2001), general equilibrium models (Gennote and Marsh 1993), and speculative storage model with rational expectations (Michealides and Ng 2000) from Economics, and compartment models from Pharmacokinetics (Mallet et al. 1988). Suppose that your data are the time series {y 1 , y 2 , ... ,y n  }, and the model that you want to estimate, or true model, is characterized by the vector of parameters θ.   The EMM procedure can be described with the following steps. Projection onto the auxiliary model. Choose an auxiliary model (or score generator) with transition density    parameterized by the pseudo parameter  , where Y t-1 = {y t -1, ... , y 1 }. The auxiliary model must approximate the true data generating process as closely as possible and be specified such that ML is feasible. The only identification requirement is that the dimension of the pseudo parameter  be greater or equal to that of the structural parameter θ. Estimate then  with the ML method. The ML estimator  satisfies the first order conditions   where   denotes the score function of the auxiliary model. Covariance matrix computation. Compute a consistent estimator of the asymptotic covariance matrix of the sample pseudo score vector. If the pseudo model is close enough to the structural model, in a suitable sense, Gallant and Long (1997) show that such an estimator may be obtained from the formula SMM. For a fixed value of the parameter vector θ, and an arbitrarily large T, simulate the series  from the true model, and find the value     that minimizes the quantity   where is the sample moment condition evaluated at the fixed estimated pseudo parameter  . Note that m T  depends on the parameter θ only through the simulated series  . The Strong Law of Large Numbers guarantees that  converges almost surely, for  , to the population moment   which, in general, cannot be computed analytically due to the presence of latent variables in the model. Under suitable conditions, the EMM estimator is consistent and asymptotically normal (Gallant and Tauchen 1996). More precisely, if   and  are the true values of the parameters of the structural and auxiliary model, respectively, then   where  A consistent estimator of the asymptotic variance-covariance matrix of   is given by     This fact can be used to compute test statistics for the EMM estimator. A Stochastic Volatility Model   Consider the following stochastic volatility model for a series of returns y t . for t = 1, ... , n, -1<b<1, s>0, and I 2  is the 2-dimensional identity matrix. The structural parameter set θ is given by (a, b, s). The data in this example consist of 1,000 observations generated from the preceding model using the following SAS statements. The values of the parameters are set to (a = -0.736, b = 0.9, s = 0.363).   /*-------------- Define n and T --------------------*/ %let nobs=1000; %let nsim=20000; /* ------------- Generate model data ---------------*/ data _tmpdata; a = -0.736; b=0.9; s=0.363; ll=sqrt( exp(a/(1-b)));; do i=-10 to &nobs; u = rannor( 101 ); z = rannor( 98761 ); lnssq = a + b*log(ll**2) +s*u; st = sqrt(exp(lnssq)); ll = st; y = st * z; if i > 0 then output; end; run;   Projection onto the GARCH(1,1) Model   Andersen et al. (1999) show that the GARCH(1,1) is an appropriate auxiliary model that leads to a good performance of the EMM estimator for such a structural model. The analytical expression of the GARCH(1,1) model with mean zero is   The pseudo parameter vector   is given by  . One advantage of such a class of models is that the conditional density ofy t is Gaussian, that is and, therefore, the score vector can easily be computed analytically. The AUTOREG procedure provides ML estimates  . The output is stored in the garchout data set, while the ML estimates are stored in the garchest data set. /*-------- estimate GARCH(1,1) model ---------------*/ proc autoreg data=_tmpdata outest=garchest; model y = / noint garch=( q=1, p=1 ) ; output out=garchout cev=gsigmasq r=resid; run;   Covariance Matrix Computation The ML estimates of the GARCH(1,1) model are used in the following SAS DATA step statements to compute the variance-covariance matrix  .   /* --------------- compute the V matrix ------------------*/ data vvalues; set garchout(keep=y gsigmasq resid); /* --- compute the scores of the GARCH model -----*/ score_1 = (-1 + y**2/gsigmasq)/ gsigmasq; score_2 = (-1 + y**2/gsigmasq)*lag(gsigmasq) / gsigmasq; score_3 = (-1 + y**2/gsigmasq)*lag(y**2) / gsigmasq; array score{*} score_1-score_3; array v_t{*} v_t_1-v_t_6; array v{*} v_1-v_6; /* - compute the external product of the score vectors - */ do i=1 to 3; do j=i to 3; v_t{j*(j-1)/2 + i} = score{i}*score{j}; end; end; /* -------- average them over t ------- */ do s=1 to 6; v{s}+ v_t{s}/&nobs; end; run;     The last observation of the V values data set contains the values of the  matrix. Obs v_1 v_2 v_3 v_4 v_5 v_6 1000 7315358.77 5282.53 4.37633 2803.39 2.36076 5.53660 Figure 1: Values of    The   matrix must be formatted to satisfy the requirements of the VDATA= option of the MODEL procedure. Please refer to the "Input Data Sets" section in "The Model Procedure" chapter of the SAS/ETS User's Guide for more information regarding the VDATA= data set. /* --------- Create a dataset acceptable to PROC MODEL --------- */ /* ---- Transpose the last obs in the dataset ----- */ proc transpose data=vvalues(firstobs=&nobs keep=v_1-v_6) out=tempv; run; /* ----- Add eq and instrument labels -----*/ data vhat; set tempv(drop=_name_); value = col1; drop col1; input _type_ $ eq_row $ eq_col $ inst_row $ inst_col $; datalines; gmm m1 m1 1 1 /* intcpt is the only instrument */ gmm m1 m2 1 1 gmm m2 m2 1 1 gmm m1 m3 1 1 gmm m2 m3 1 1 gmm m3 m3 1 1 ; run;   The vhat data set contains the formatted matrix. Obs value _type_ eq_row eq_col inst_row inst_col 1 7315358.77 gmm m1 m1 1 1 2 5282.53 gmm m1 m2 1 1 3 4.38 gmm m2 m2 1 1 4 2803.39 gmm m1 m3 1 1 5 2.36 gmm m2 m3 1 1 6 5.54 gmm m3 m3 1 1 Figure 2: Formatted    SMM The following statements generate two independent random sequences of length T=20,000 for the errors and utilize them to find the EMM estimates for the structural model using the GMM option of the FIT statement in the MODEL procedure. Please note that the   matrix serves as weighting matrix in the VDATA= option, and that the scores of the GARCH(1,1) auxiliary model evaluated at the ML estimates are the moment conditions in the GMM step. Since the number of parameters to estimate (3) is equal to the number of moment equations (3) times the number of instruments (1), the model is exactly identified and the objective function will have value zero at the minimum. /* ---------------------- FIND EMM ESTIMATES ---------------------*/ /* ------- Simulate the error sequences with datastep -------*/ data emm; set garchest(obs=1 keep = _mse_ _ah_0 _ah_1 _gh_1); do i=1 to &nsim; u = rannor( 8801 ); z = rannor( 9701 ); output; end; drop i; run; /* ------------- Estimate parameters with GMM -----------------*/ title 'EMM Estimates'; proc model data=emm; parms a -.7 b .9 s .363; instrument _exog_ / intonly; /*-- Describe the structural model ----*/ lsigmasq = xlag(sigmasq,exp(a)); lnsigmasq = a + b * log(lsigmasq) + s * u; sigmasq = exp( lnsigmasq ); ysim = sqrt(sigmasq) * z; /*--- Volatility of the GARCH model ---*/ gsigmasq = _ah_0 + _gh_1*xlag(gsigmasq, _mse_) + _ah_1*xlag(ysim**2, _mse_); /*---- Scores of the GARCH model as moment conditions ----*/ eq.m1 = (-1 + ysim**2/gsigmasq)/ gsigmasq; eq.m2 = (-1 + ysim**2/gsigmasq)*xlag(gsigmasq, _mse_) / gsigmasq; eq.m3 = (-1 + ysim**2/gsigmasq)*xlag(ysim**2, _mse_) / gsigmasq; /*----- Fit the GARCH scores using GMM and estimated Vhat -----*/ fit m1 m2 m3 / gmm vdata=vhat kernel=(bart,0,) no2sls; bounds s > 0, -1<b<1; run;   The output of the MODEL procedure is as follows:   EMM Estimates The MODEL Procedure Model Summary Parameters 3 Equations 3 Number of Statements 12 Parameters(Value) a(-0.7) b(0.9) s(0.363) Equations m1 m2 m3 The 3 Equations to Estimate m1 = F(a, b, s) m2 = F(a, b, s) m3 = F(a, b, s) Instruments 1   EMM Estimates The MODEL Procedure Nonlinear GMM Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| a -0.50466 0.0778 -6.49 <.0001 b 0.928844 0.0107 87.05 <.0001 s 0.294869 0.0304 9.71 <.0001 Figure 3: PROC MODEL Output Computational Considerations The SMM portion of the EMM estimation procedure relies on the simulated series  . Different random seeds will generate different series, and, therefore, for any finite T, the EMM estimates will vary according to the random seed chosen. However, for T sufficiently large, this variation is negligible when compared to the variation due to the randomness in the sample. For example, using the code provided in Appendix A, you can see that the standard errors of the EMM estimates for T=20,000, computed with a Monte Carlo experiment using 500 different random seeds for the simulated series  , are 0.0052 for a, 0.00071 for b, and 0.0023 for s. You can compare these values to the standard errors due to random variations in the sample as computed in Andersen et al. (1999): 0.60 for a, 0.88 for b, and 0.20 for s. As with all nonlinear estimation problems, the computational time required by the EMM estimation can vary a great deal depending on several factors, particularly the choice of the initial values. If the initial values are sufficiently close to the solution, it may take as few as 30 CPU seconds on a Pentium PIII 500Mhz class computer. If the starting values are far from the solution, it may take a several minutes, or not converge at all, depending on the starting values themselves and on the optimization options in PROC MODEL. For computational convenience, the starting values in the above example are the "true" values of the parameters. It is, however, advisable to choose a grid of starting values and repeat the estimation for each set in the grid. Unfortunately, this task is computationally burdensome, and cannot be accomplished using the START= and STARTITER= options of the FIT statement in PROC MODEL, because such options use the 2SLS method and not the moment conditions provided. If you have access to several computers running SAS, you may use the %Distribute macro of Doninger and Tobias (2001) to accomplish such a task in parallel on all the computers and reduce significantly the overall computational time. The code in Appendix A can be easily modified to this purpose.   Appendix A   This appendix contains the code for a Monte Carlo study on the properties of the EMM estimator for the stochastic volatility model. The series     is simulated using 500 different seeds for each length T, starting from T=20,000 to T=200,000, at intervals of 20,000, and the EMM estimates are computed for each generated series. The distributions of the resulting estimates are summarized by series length using PROC MEANS. This task is accomplished with parallel computing on seven PCs thanks to the %Distribute macro described in Doninger and Tobias (2001). User names, passwords, and computer addresses have been altered to protect confidentiality.   /***************************************************************/ /* Filename: EMM_DISTR.SAS */ /* Purpose: study distribution of EMM estimates using parallel */ /* distributed computing with %Distribute */ /* Date: 10/03/01 */ /* Modified: */ /* Author: Snow White */ /***************************************************************/ /* / Define libraries where to save results. /------------------------------------------------------------------*/ libname EmmDist "u:\distrib"; /* for pcs */ /* / Include %Distribute and other relevant macros / -----------------------------------------------------------------*/ %inc "u:\distrib\Distribute.sas" / nosource; /* --------------------- UNDISTRIBUTED PART ---------------------- */ /* --------- MACRO for DATA, GARCH ESTIMATES AND VHAT ---------- */ %macro gen_data( nobs = 1000, /* number of obs in the observed data set */ vhat = EmmDist.vhat, /* name of the data set for Vhat */ garchest = EmmDist.garchest, /* data set for GARCH estimates */ seed = 101 /* seed for u and z */ ); data _tmpdata; a = -0.736; b=0.9; s=0.363; ll=sqrt( exp(a/(1-b)));; do i=-10 to &nobs; u = rannor( &seed ); z = rannor( &seed ); lnssq = a + b*log(ll**2) +s*u; st = sqrt(exp(lnssq)); ll = st; y = st * z; if i > 0 then output; end; run; /*-------- estimate GARCH(1,1) model ---------------*/ proc autoreg data=_tmpdata outest=&garchest; model y = / noint garch=( q=1, p=1 ) ; output out=garchout cev=gsigmasq r=resid; run; quit; /* --------------- compute the V matrix ------------------*/ data vvalues; set garchout(keep=y gsigmasq resid); /* --- compute the scores of the GARCH model -----*/ score_1 = (-1 + y**2/gsigmasq)/ gsigmasq; score_2 = (-1 + y**2/gsigmasq)*lag(gsigmasq) / gsigmasq; score_3 = (-1 + y**2/gsigmasq)*lag(y**2) / gsigmasq; array score{*} score_1-score_3; array v_t{*} v_t_1-v_t_6; array v{*} v_1-v_6; /* --- compute the external product of the score vectors -- */ do i=1 to 3; do j=i to 3; v_t{3*(i>1) + 2*(i>2)+j-i+1} = score{i}*score{j}; end; end; /* -------- average them over t ------- */ do s=1 to 6; v{s}+ v_t{s}/&nobs; end; run; quit; /* ---------- Create a data set acceptable to PROC MODEL -------- */ /* ---- Transpose the last obs in the data set ----- */ proc transpose data=vvalues(firstobs=&nobs keep=v_1-v_6) out=tempv; run; quit; /* ----- Add eq and instrument labels -----*/ data labels; _type_ = 'gmm'; inst_row ='1'; inst_col ='1'; do i=1 to 3; do j=i to 3; eq_row=compress("m"||j); eq_col=compress("m"||i); output; end; end; drop i j; run; data &vhat; merge labels tempv(drop=_name_); value=col1; drop col1; run; %mend; /* ------------------- END UNDISTRIBUTED PART ------------------- */ /* / Create pc host data set for pcs /------------------------------------------------------------------*/ data _hosts; script= "c:\Program Files\SAS Institute\sas\connect\saslink\tcpwin.scr"; input host $; datalines; Doc Sleepy Sneezy Bashful Happy Grumpy Dopey ; run; /* / Login to hosts. /------------------------------------------------------------------*/ %let USER = SnowWhite; %let PASS = apple; %SignOn; /* / %RINIT macro /-------------------------------------------------------------- */ %macro RInit; rsubmit remote=Host&iHost wait=yes macvar=Status&iHost; %macro FirstRSub; *options nonotes; /* create empty output data sets */ data REmm_est; if (0); run; data REmm_ods; if (0); run; /* -------- Rsubmit the MACRO for EMM ESTIMATION ------- */ %macro Emm( garchest=garchest, /* GARCH estimates */ vhat = vhat, /* data set of Vhat */ Emm_est = REmm_est, /* output data set of estimates */ Emm_ods = REmm_ods, /* ODS ParameterEstimates data set */ seed = &rem_Seed, /* seed for u and z */ nsim = &rem_Nsim /* number of simulated obs */ ); /* -- Simulate error sequences with datastep --*/ data _Emm; set &garchest(obs=1 keep = _mse_ _ah_0 _ah_1 _gh_1); do i=1 to &nsim; u = rannor( &seed ); z = rannor( &seed ); output; end; drop i; run; /* ------ Estimate parameters with GMM ---------*/ ods output ParameterEstimates=_Emm_ods; ods output show; /* save the par estimates with ods */ proc model data=_Emm outparms=_Emm_est; parms a -.7 b .9 s .363; instrument _exog_ / intonly; /*-- Describe the structural model ----*/ lsigmasq = xlag(sigmasq,exp(a)); lnsigmasq = a + b * log(lsigmasq) + s * u; sigmasq = exp( lnsigmasq ); ysim = sqrt(sigmasq) * z; /*--- Volatility of the GARCH model ---*/ gsigmasq = _ah_0 + _gh_1*xlag(gsigmasq, _mse_) + _ah_1*xlag(ysim**2, _mse_); /*- Scores of the GARCH model as moment conditions -*/ eq.m1 = (-1 + ysim**2/gsigmasq)/ gsigmasq; eq.m2 = (-1 + ysim**2/gsigmasq)*xlag(gsigmasq,_mse_) / gsigmasq; eq.m3 = (-1 + ysim**2/gsigmasq)*xlag(ysim**2,_mse_) / gsigmasq; /*- Fit the GARCH scores using GMM and estimated Vhat -*/ fit m1 m2 m3 / gmm vdata=vhat kernel=(bart,0,); bounds s > 0, -1<b<1; run; ods output close; /* -- add relevant info in Estimates data set -- */ data &Emm_est; set &Emm_est _Emm_est; nsim = &nsim; seed = &seed; run; /* -- add relevant info in ODS output data set -- */ data &Emm_ods; set &Emm_ods _Emm_ods; nsim = &nsim; seed = &seed; run; %mend; %mend; %macro TaskRSub; %Emm; /* Run 1 EMM estimation at each remote submission */ %mend; endrsubmit; %mend; /* / Loop over different simulated series lengths /------------------------------------------------------------------*/ %macro Emm_Sim(start=20000, end=200000, by=20000 ); data EmmDist.Emm_est; if (0); run; /* create empty data sets */ data EmmDist.Emm_ods; if (0); run; /* for output */ %do nsim=&start %to &end %by &by; %Distribute; %RCollect(REmm_est, REmm_est); data EmmDist.Emm_est; set EmmDist.Emm_est REmm_est; run; %RCollect(REmm_ods, REmm_ods); data EmmDist.Emm_ods; set EmmDist.Emm_ods REmm_ods; run; %end; %mend; /* / Define remote length of simulated serie. /------------------------------------------------------------------*/ %macro Macros; %syslput rem_Nsim=&nsim; %mend; /* / Generate observed data and transfer to remote hosts Work dir /------------------------------------------------------------------*/ %gen_data; %macro trandata; %do i=1 %to 7; data RWork%left(&i).garchest; set EmmDist.garchest; run; data RWork%left(&i).vhat; set EmmDist.vhat; run; %end; %mend; %trandata; /* / Arguments for the distribution. /------------------------------------------------------------------*/ %let NIter = 1; %let NIterAll = 500; %let Debug = 0; /* / Run all the stuff and analyze results. /------------------------------------------------------------------*/ %Emm_Sim(start=20000, end=200000, by=20000); proc means data=EmmDist.Emm_est n mean stderr lclm uclm; var a b s; by nsim; run; %SignOff;     References   Andersen, T.G., Chung, H-J. and Sorensen, B.E. (1999), "Efficient Method of Moments Estimation of a Stochastic Volatility Model: A Monte Carlo Study," Journal of Econometrics, 91, 61-87. Andersen, T.G. and Sorensen, B.E. (1996), "GMM Estimation of a Stochastic Volatility Model: A Monte Carlo Study," Journal of Business and Economic Statistics, 14, 328-352. Bansal, R., Gallant, A.R., Hussey, R., Tauchen, G.E. (1993), "Computational Aspects of Nonparametric Simulation Estimation." In Belsey, D.A. (Ed.), Computational Techniques for Econometrics and Economic Analysis. Boston, MA: Kluwer Academic Publishers, 3-22. Bansal, R., Gallant, A.R., Hussey, R., Tauchen, G.E. (1995), "Nonparametric Estimation of Structural Models for High-Frequency Currency Market Data," Journal of Econometrics, 66, 251-287. Doninger, C. and Tobias, R. (2001), "The %Distribute System for Large-Scale Parallel Computation in the SAS System." [http://www.sas.com/rnd/app/papers/distConnect.pdf] . Cary, NC: SAS Institute, Inc. Accessed 12 October 2001. Duffie, D. and Singleton, K.J. (1993), "Simulated Moments Estimation of Markov Models of Asset Prices," Econometrica 61, 929-952. Gallant, A.R. and Long, J. (1997), "Estimating Stochastic Differential Equations Efficiently by Minimum Chi-squared," Biometrika, 84, 125-141. Gallant, A.R. and Tauchen, G.E. (1996), "Which Moments to Match?" Econometric Theory, 12, 657-681. Gallant, A.R. and Tauchen, G.E. (2001), "Efficient Method of Moments," Working Paper. [http://www.econ.duke.edu/ get/wpapers/ee.pdf] , accessed 12 September 2001. Gennotte, G. and Marsh, T.A. (1993), "Variations in Economic Uncertainty and Risk Premiums on Capital Assets," European Economic Revue, 37, 1021-1041. Mallet, A., Mentré, F., Steimer, J-L., and Lokiec, F. (1988), "Nonparametric Maximum Likelihood Estimation for Population Pharmacokinetics, with Application to Cyclosporine," Journal of Pharmacokinetics and Biofarmaceutics, 16, 311-327. Michaelides, A. and Ng, S. (2000), "Estimating the Rational Expectations Model of Speculative Storage: A Monte Carlo Comparison of Three Simulation Estimators," Journal of Econometrics, 96(2), 231-266. Olsen, L.F. and Schaffer, W.M. (1990), "Chaos versus Noisy Periodicity: Alternative Hypotheses for Childhood Epidemics," Science, 249, 499-504. SAS Institute, Inc. (1999), SAS/ETS User's Guide, Version 8, Volumes 1 and 2, Cary, NC: SAS Institute, Inc. ... View more

Overview   The transcendental logarithmic (translog) functional form can capture many of the attributes of a cost function that are implied by economic theory. Because of this flexibility, it has been widely used for studying production relationships. Before logarithms are taken, the function is The function takes w and y as its inputs and returns minimum cost, where w is a vector of prices for the inputs to production and y is a single output. N is the total number of inputs and the α's are the parameters of the function. It may not be apparent from the equation above, but the translog form nests the Cobb-Douglas form and is equivalent if all α i j  and α iy  are zero. One disadvantage of the previous form is that it is not linear in parameters. A standard technique when dealing with power functions like the translog cost function is to take logarithms. The resulting function is linear in parameters and standard statistical techniques can be used for estimation. After taking logarithms the function is   While it is possible to include terms to account for technological progress, the specification used here assumes that cost is independent of time. Using Shepherd’s lemma, the derived demand equations are where   is the cost share of the i th input. The cost function is assumed to be continuous, so Young’s Theorem concerning symmetry of the second derivatives restricts α i j   = α ji   for all i ≠j.   The result of this derivation is a system of N + 1 equations consisting of N derived demand equations and one cost function. Homogeneity of the first degree implies for all i and j.   It is also possible to impose constant returns to scale – equivalent to imposing homogeneity in y – and details of this procedure can be found in Diewert and Wales (1987). Global concavity can also be imposed on this specification by forcing the matrix [α i j ] to be negative semidefinite. A technique for accomplishing this can be found in Jorgenson (1986). This procedure tends to restrict the price elasticities in an undesirable way and destroys the flexibility of the form.   Analysis   This example shows the estimation of a system of nonlinear derived demand equations based on a translog cost function using iterated seemingly unrelated regression (ITSUR). Data pertains to the United States textiles manufacturing sector (standard industrial classification code 22). The series runs from 1949 to 2001 and contains real quantity indices, real price indices, and real cost measures for a single aggregate industry output and five aggregate inputs: capital (K), labor (L), energy (E), materials (M), and services (S). The original dataset, and information on other industrial sectors, can be obtained from the Multifactor Productivity Home Page of the Bureau of Labor Statistics at http://www.bls.gov/mfp/. Derived demand equations are each expressed with a cost share as the endogenous variable. Because this dataset does not contain explicit information on cost shares, the shares must be formed by taking the ratio of the value of each input and the cost measure.   data klems; set klems; array values {5} vk vl ve vm vs; array costshares {5} sk sl se sm ss; cost = sum(vk,vl,ve,vm,vs); do i = 1 to 5; costshares{i} = values{i}/cost; end; run; A plot of the cost shares over time is produced by the following statements.   proc sgplot data = klems; series x = year y = sk / markers markerattrs =(symbol=circle); series x = year y = sl / markers markerattrs =(symbol=square); series x = year y = se / markers markerattrs =(symbol=star); series x = year y = sm / markers markerattrs =(symbol=diamond); series x = year y = ss / markers markerattrs =(symbol=hash); title 'Factor Cost Shares'; yaxis label = 'Cost Share'; run;   Figure 1: Change in Cost Shares   One benefit of the translog form is that the system of factor demands gives nearly the same information as the cost function. The only parameter of the cost function that is not captured by the derived demand system is the intercept term which is not used in calculating price or substitution elasticities. To be as parsimonious as possible, the cost function is typically not estimated. In this example, the MODEL procedure is used to fit four derived demand equations. One of the equations (the derived demand equation for services) has been arbitrarily dropped from estimation; only N - 1 of the factor demands are linearly independent because the dependent variables are cost shares which must sum to one. The parameters of the dropped derived demand equation can be recovered after estimation through homogeneity and symmetry restrictions. The following statements estimate the system of derived demand equations without imposing any restrictions from theory. Likelihood ratio tests are conducted to determine whether homogeneity and symmetry are satisfied both singularly and jointly.   proc model data = klems; parameters a_k gkk gkl gke gkm gks gky a_l glk gll gle glm gls gly a_e gek gel gee gem ges gey a_m gmk gml gme gmm gms gmy; endogenous sk sl se sm; exogenous pk pl pe pm ps y; /*System of Derived Demand Equations*/ sk = a_k + gkk*log(pk) + gkl*log(pl) + gke*log(pe) + gkm*log(pm) + gks*log(ps) + gky*log(y); sl = a_l + glk*log(pk) + gll*log(pl) + gle*log(pe) + glm*log(pm) + gls*log(ps) + gly*log(y); se = a_e + gek*log(pk) + gel*log(pl) + gee*log(pe) + gem*log(pm) + ges*log(ps) + gey*log(y); sm = a_m + gmk*log(pk) + gml*log(pl) + gme*log(pe) + gmm*log(pm) + gms*log(ps) + gmy*log(y); fit sk sl se sm / itsur; test "Homogeneity" gkk+gkl+gke+gkm+gks=0, glk+gll+gle+glm+gls=0, gek+gel+gee+gem+ges=0, gmk+gml+gme+gmm+gms=0, / lr; test "Symmetry" gkl=glk, gke=gek, gkm=gmk, glm=gml, gle=gel, gem=gme, / lr; test "Joint Homogeneity and Symmetry" gkk+gkl+gke+gkm+gks=0, glk+gll+gle+glm+gls=0, gek+gel+gee+gem+ges=0, gmk+gml+gme+gmm+gms=0, gkl=glk, gke=gek, gkm=gmk, glm=gml, gle=gel, gem=gme, / lr; run;   Because the form of the elasticities is somewhat complicated, it can be difficult to interpret the values and signs of parameter estimates. The test results in Figure 2  show that both symmetry and homogeneity are rejected. A common practice is to assume that such restrictions hold and to impose them in estimation.   Figure 2: Tests of Symmetry and Homogeneity   The MODEL Procedure   Test Results Test Type Statistic Pr > ChiSq Label Homogeneity L.R. 114.22 <.0001 gkk+gkl+gke+gkm+gks=0, glk+gll+gle+glm+gls=0, gek+gel+gee+gem+ges=0, gmk+gml+gme+gmm+gms=0 Symmetry L.R. 109.72 <.0001 gkl=glk, gke=gek, gkm=gmk, glm=gml, gle=gel, gem=gme Joint Homogeneity and Symmetry L.R. 240.80 <.0001 gkk+gkl+gke+gkm+gks=0, glk+gll+gle+glm+gls=0, gek+gel+gee+gem+ges=0, gmk+gml+gme+gmm+gms=0, gkl=glk, gke=gek, gkm=gmk, glm=gml, gle=gel, gem=gme   The following code imposes both symmetry and homogeneity restrictions on the underlying model. A likelihood ratio test of constant returns to scale is introduced and a Chow test is added to the FIT statement. proc model data = klems; parameters a_k gkk gkl gke gkm gks gky a_l glk gll gle glm gls gly a_e gek gel gee gem ges gey a_m gmk gml gme gmm gms gmy; endogenous sk sl se sm; exogenous pk pl pe pm ps y; restrict /*Homogeneity Restrictions*/ gks=0-gkk-gkl-gke-gkm, gls=0-gkl-gll-gle-glm, ges=0-gke-gle-gee-gem, gms=0-gkm-glm-gem-gmm, /*Symmetry Restrictions*/ gkl=glk, gke=gek, gkm=gmk, gle=gel, glm=gml, gem=gme; /*System of Derived Demand Equations*/ sk = a_k + gkk*log(pk) + gkl*log(pl) + gke*log(pe) + gkm*log(pm) + gks*log(ps) + gky*log(y); sl = a_l + glk*log(pk) + gll*log(pl) + gle*log(pe) + glm*log(pm) + gls*log(ps) + gly*log(y); se = a_e + gek*log(pk) + gel*log(pl) + gee*log(pe) + gem*log(pm) + ges*log(ps) + gey*log(y); sm = a_m + gmk*log(pk) + gml*log(pl) + gme*log(pe) + gmm*log(pm) + gms*log(ps) + gmy*log(y); fit sk sl se sm / itsur chow = (24) outest=est; test "Constant Returns to Scale" gky=0, gly=0, gey=0, gmy=0, / lr; run; The symmetry restriction shrinks the number of parameters of the model considerably. This is particularly useful in cases where the time series is not long and degrees of freedom need to be conserved. Figure 3 shows the parameter estimates of the MODEL procedure.   Figure 3: The Restricted Model The MODEL Procedure   Nonlinear ITSUR Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| Label a_k 0.109996 0.0219 5.02 <.0001 SK Intercept gkk 0.062014 0.00357 17.38 <.0001 SK K Price gkl -0.01898 0.00725 -2.62 0.0118 SK L Price gke -0.00179 0.000836 -2.15 0.0369 SK E Price gkm -0.03872 0.00610 -6.35 <.0001 SK M Price gks -0.00252 0.00342 -0.74 0.4648 SK S Price gky -0.00104 0.00496 -0.21 0.8354 SK Output a_l 0.865473 0.0903 9.58 <.0001 SL Intercept glk -0.01898 0.00725 -2.62 0.0118 SL K Price gll -0.00999 0.0360 -0.28 0.7826 SL L Price gle -0.00106 0.00420 -0.25 0.8013 SL E Price glm -0.06766 0.0248 -2.73 0.0087 SL M Price gls 0.097692 0.0214 4.57 <.0001 SL S Price gly -0.11488 0.0200 -5.74 <.0001 SL Output a_e 0.012747 0.00984 1.30 0.2013 SE Intercept gek -0.00179 0.000836 -2.15 0.0370 SE K Price gel -0.00106 0.00420 -0.25 0.8013 SE L Price gee 0.029876 0.00112 26.76 <.0001 SE E Price gem -0.01954 0.00275 -7.12 <.0001 SE M Price ges -0.00748 0.00325 -2.30 0.0257 SE S Price gey 0.005808 0.00217 2.68 0.0101 SE Output a_m -0.08173 0.0701 -1.17 0.2492 SM Intercept gmk -0.03872 0.00610 -6.35 <.0001 SM K Price gml -0.06766 0.0248 -2.73 0.0087 SM L Price gme -0.01954 0.00275 -7.12 <.0001 SM E Price gmm 0.146849 0.0222 6.62 <.0001 SM M Price gms -0.02092 0.0103 -2.03 0.0477 SM S Price gmy 0.113729 0.0157 7.27 <.0001 SM Output Restrict0 -563.453 242.8 -2.32 0.0187 gks=0-gkk-gkl-gke-gkm Restrict1 82.23307 193.0 0.43 0.6747 gls=0-gkl-gll-gle-glm Restrict2 321.7446 689.2 0.47 0.6455 ges=0-gke-gle-gee-gem Restrict3 -279.71 211.0 -1.33 0.1879 gms=0-gkm-glm-gem-gmm Restrict4 -261.228 196.3 -1.33 0.1860 gkl=glk Restrict5 -1041.84 755.0 -1.38 0.1700 gke=gek Restrict6 -27.855 219.3 -0.13 0.9005 gkm=gmk Restrict7 -880.821 742.7 -1.19 0.2396 gle=gel Restrict8 259.513 215.6 1.20 0.2326 glm=gml Restrict9 1103.343 320.8 3.44 0.0003 gem=gme   The majority of the parameter estimates are significant. Insignificant parameters are statistically equivalent to zero and imply that corresponding elasticities of substitution are equal to the Cobb-Douglas value of one. The TEST statement tests whether this industry exhibits constant returns to scale in the range of the sample. The CHOW option in the FIT statement performs a Chow test for a structural break at the twenty-fourth year of the sample. The twenty-fourth year of the sample is 1973; in October of that year OPEC declared an oil embargo. Markets were affected by significant shocks to oil prices and gasoline was rationed in the United States. The results of the two tests can be seen in Figure 4.   Figure 4: CRS and Chow Test Results Test Results Test Type Statistic Pr > ChiSq Label Constant Returns to Scale L.R. 74.27 <.0001 gky=0, gly=0, gey=0, gmy=0 Structural Change Test Test Break Point Num DF Den DF F Value Pr > F Chow 24 23 2 0.28 0.9560   The null hypothesis of constant returns to scale is rejected. The null hypothesis of the Chow Test cannot be rejected. Even with the turmoil of the oil embargo, there is no evidence of a structural break in 1973. Parameter estimates are stored in the est dataset and used in calculating price elasticities and elasticities of substitution in the example "Calculating Elasticities from a Translog Cost Function." References   Blackorby, C., and Russell, R. R. (1989). “Will the Real Elasticity of Substitution Please Stand Up? (A Comparison of the Allen/Uzawa and Morishima Elasticities).” American Economic Review 79:882–888. Chambers, R. G. (1988). Applied Production Analysis: A Dual Approach. New York: Cambridge University Press. Diewert, W. E., and Wales, T. J. (1987). “Flexible Functional Forms and Global Curvature Conditions.” Econometrica 55:43–68. Jorgenson, D. (1986). “Econometric Methods for Modeling Producer Behavior.” In Handbook of Econometrics, edited by Z. Griliches, and M. D. Intriligator, 1841–1915. Amsterdam: North-Holland. ... View more

Overview   The Almost Ideal Demand System (AIDS) model of Deaton and Muellbauer (1980b) has enjoyed great popularity in applied demand analysis. Starting from a specific cost function, the AIDS model gives the share equations in an n-good system as   where w i is the share associated with the i th good,  is the constant coefficient in the i th share equation,  is the slope coefficient associated with the j th good in the i th share equation, pj is the price on the j th good. X is the total expenditure on the system of goods given by in which q i is the quantity demanded for the i th good. P is the price index defined by in the nonlinear AIDS model. Deaton and Muellbauer (1980a) also suggested a linear approximation of the nonlinear AIDS model by specifying a linear price index given by that gives rise to the linear approximate AIDS (LA-AIDS) model. In practice, the LA-AIDS model is more frequently estimated than the nonlinear AIDS model. Conservation implies the following restrictions on the parameters in the nonlinear AIDS model: Homogeneity is satisfied if and only if, for all i  and symmetry is satisfied if  One advantage of the AIDS model is that the homogeneity and symmetry restrictions are easily imposed and tested. Analysis   In this example, you can estimate a nonlinear AIDS model and LA-AIDS model using U.S. meat products data. The data consists of aggregate quarterly retail price and per capita consumption for four meat product categories including beef, pork, chicken, and turkey. The data period covers the first quarter of 1974 to the last quarter of 1999. The data were obtained from various USDA sources listed in the references. More detailed description of the data construction may be found in Piggott (1997). To incorporate seasonality and trend in the U.S. meat consumption data, you augment the AIDS model with trigonometric variables and a time trend variable. Thus the share equations to be estimated are expressed as follows:  where aci and asi represent parameters on the trigonometric variables, and ati is the parameter on the time trend variable. Since the budget shares sum up to 1, these parameters should satisfy   The following is a portion of the code used to create the data set. A second set of SAS statements is used to create the time trend variable and the seasonal variables. Total expenditures, X, and shares for each good are also computed.     data aids_; input year qtr pop beef_q pork_q chick_q turkey_q beef_p pork_p chick_p turkey_p cpi pc_exp; date = yyq(year,qtr); format date yyq6.; ... more datalines ... ; run;   data aids; set aids_; if year < 1975 then delete; t = _n_ ; co1 = cos(1/2*3.14159*t); co2 = cos(2/2*3.14159*t); si1 = sin(1/2*3.14159*t); si2 = sin(2/2*3.14159*t); x = beef_p*beef_q + pork_p*pork_q + chick_p*chick_q + turkey_p*turkey_q; w_beef = beef_p*beef_q/x; w_pork = pork_p*pork_q/x; w_chick = chick_p*chick_q/x; w_turkey = turkey_p*turkey_q/x; run;   The plot of the budget shares found in Figure 1 .   proc gplot data=aids; plot w_beef*date w_pork*date w_chick*date w_turkey*date/ overlay cframe=ligr haxis=axis1 vaxis=axis2; title 'Budget Shares Plots'; footnote1 c=blue ' * beef ' c=red ' . pork ' c=green ' o chicken ' c=black ' + turkey '; symbol1 c=blue i=join v=star; symbol2 c=red i=join v=dot; symbol3 c=green i = join v = circle ; symbol4 c=black i = join v = plus ; axis1 label=('Time') ; axis2 label=(angle=90 'Budget Share'); run; quit;   Figure 1: Budget Shares Plots Price and quantity variables in the data set are labeled in an obvious manner. Three variable names in the DATA step require some explanation. The variable named pop contains the U.S. population, cpi contains the consumer price index for all goods, and pc_exp is the per capita nominal personal consumption expenditures on all goods. To obtain mean-scaled price data, you use the MEANS procedure. The means of each of the price variables are written to the OUT= data set as the b_m, p_m, c_m, and t_m variables. The original price variables are then divided by the mean variables in a subsequent DATA step.   data aids ; if _n_ = 1 then set pm ; set aids ; pb_ = (beef_p/b_m); pp_ = (pork_p/p_m); pc_ = (chick_p/c_m); pt_ = (turkey_p/t_m); lpb = log(pb_); lpp = log(pp_); lpc = log(pc_); lpt = log(pt_); lx = log(x); /* Specify the linear price index for the LA-AIDS model */ lp0 = sum(w_beef*lpb + w_pork*lpp + w_chick*lpc + w_turkey*lpt); lxp = lx-lp0; run;   The LA-AIDS model is estimated using the SYSLIN procedure. The ITSUR option specifies the iterated seemingly unrelated regression as the estimation method. The equations to be fit are then specified. Since budget shares sum to 1 in the system, one of the share equations is deleted to deal with the singularity problem. In this example, you eliminate the share equation for turkey. Whichever one is eliminated should not have any effect on the results. The parameters associated with the share equation that is deleted can be recovered through the parameter restrictions implied by the homogeneity, symmetry, and conservation properties. The parameter restrictions can be imposed through the SRESTRICT statement.     proc syslin itsur data = aids outest = fin1; b: model w_beef = lpb lpp lpc lpt lxp co1 si1 t ; p: model w_pork = lpb lpp lpc lpt lxp co1 si1 t ; c: model w_chick = lpb lpp lpc lpt lxp co1 si1 t ; srestrict /*symmetry restrictions */ b.lpp = p.lpb , b.lpc = c.lpb , p.lpc = c.lpp , /* homogeneity restrictions*/ b.lpb + b.lpp + b.lpc + b.lpt = 0 , p.lpb + p.lpp + p.lpc + p.lpt = 0 , c.lpb + c.lpp + c.lpc + c.lpt = 0; run;   The parameter estimates in the beef equation are given in Figure 2.   Parameter Estimates - Beef Variable DF Estimate StdErr tValue Probt Intercept 1 0.397899 0.400708 0.99 0.3234 lpb 1 0.041175 0.012686 3.25 0.0016 lpp 1 -0.00407 0.008127 -0.50 0.6181 lpc 1 -0.04926 0.00541 -9.10 <.0001 lpt 1 0.012151 0.009613 1.26 0.2095 lxp 1 0.023424 0.044316 0.53 0.5984 co1 1 -0.01631 0.001631 -10.00 <.0001 si1 1 -0.00228 0.001571 -1.45 0.1497 t 1 -0.00148 0.000049 -30.45 <.0001 Figure 2: Parameter Estimates for Beef The parameter estimates in the pork equation are given in Figure 3. Parameter Estimates - Pork Variable DF Estimate StdErr tValue Probt Intercept 1 0.240957 0.28773 0.84 0.4046 lpb 1 -0.00407 0.008127 -0.50 0.6181 lpp 1 0.044733 0.008811 5.08 <.0001 lpc 1 -0.03074 0.006022 -5.10 <.0001 lpt 1 -0.00993 0.007543 -1.32 0.1915 lxp 1 0.003386 0.031819 0.11 0.9155 co1 1 0.009611 0.001187 8.10 <.0001 si1 1 0.006454 0.001137 5.68 <.0001 t 1 0.000278 0.000036 7.70 <.0001 Figure 3: Parameter Estimates for Pork The parameter estimates in the chicken equation are given in Figure 4. Parameter Estimates - Chicken Variable DF Estimate StdErr tValue Probt Intercept 1 0.579604 0.17659 3.28 0.0015 lpb 1 -0.04926 0.00541 -9.10 <.0001 lpp 1 -0.03074 0.006022 -5.10 <.0001 lpc 1 0.101387 0.007607 13.33 <.0001 lpt 1 -0.02139 0.006996 -3.06 0.0029 lxp 1 -0.05355 0.019548 -2.74 0.0074 co1 1 -0.00567 0.0007 -8.10 <.0001 si1 1 -0.00084 0.000664 -1.26 0.2114 t 1 0.000893 0.000026 33.86 <.0001 Figure 4: Parameter Estimates for Chicken The parameter restrictions estimates are reported in Figure 5: Parameter Restrictions Estimates Variable DF Estimate StdErr tValue Probt RESTRICT -1 -233.279 65.95364 -3.54 0.0003 RESTRICT -1 -130.174 75.45235 -1.73 0.0845 RESTRICT -1 -14.1833 95.49745 -0.15 0.8829 RESTRICT -1 158.8543 73.2354 2.17 0.0293 RESTRICT -1 -36.2377 88.88858 -0.41 0.6859 RESTRICT -1 -301.062 110.8651 -2.72 0.0059 Figure 5: Parameter Restrictions Estimates For the nonlinear AIDS model, you estimate the model using the MODEL procedure. Homogeneity restrictions can easily be imposed using the RESTRICT statement. The symmetry restrictions can be imposed by using the same parameter name for the corresponding variables. The nonlinear price index needs to be specified as part of the model program. The FIT statement estimates model parameters by fitting the specified equations to the data. The ITSUR option in the FIT statement specifies the estimation method as the iterated seemingly unrelated regression method. The parameters to be estimated are specified in the PARMS statement.   /* Full Nonlinear AIDS Model */ proc model data=aids; /* imposing homogeneity and symmetry and adding-up restrictions */ restrict gbb + gbp + gbc + gbt = 0 , gbp + gpp + gpc + gpt = 0 , gbc + gpc + gcc + gct = 0 , gbt + gpt + gct + gtt = 0 , ab + ap + ac + at = 1 ; /*non-linear price index*/ a0 = 0; /* restrict the constant term in the nonlinear price index to be zero */ p = a0 + ab*lpb + ap*lpp + ac*lpc + at*lpt + .5*(gbb*lpb*lpb + gbp*lpb*lpp + gbc*lpb*lpc + gbt*lpb*lpt + gbp*lpp*lpb + gpp*lpp*lpp + gpc*lpp*lpc + gpt*lpp*lpt + gbc*lpc*lpb + gpc*lpc*lpp + gcc*lpc*lpc + gct*lpc*lpt + gbt*lpt*lpb + gpt*lpt*lpp + gct*lpt*lpc + gtt*lpt*lpt ); /*share equation*/ w_beef = ab + gbb*lpb + gbp*lpp + gbc*lpc + gbt*lpt + bb*(lx-p) + abco1*co1 + absi1*si1 + ab_t*t ; w_pork = ap + gbp*lpb + gpp*lpp + gpc*lpc + gpt*lpt + bp*(lx-p) + apco1*co1 + apsi1*si1 + ap_t*t ; w_chick = ac + gbc*lpb + gpc*lpp + gcc*lpc + gct*lpt + bc*(lx-p) + acco1*co1 + acsi1*si1 + ac_t*t ; fit w_beef w_pork w_chick / itsur nestit outs=rest estdata=fin0 outest=fin2 out = resid2 converge = .00001 maxit = 1000 ; parms ab bb gbb gbp gbc gbt abco1 absi1 ab_t ap bp gpp gpc gpt apco1 apsi1 ap_t ac bc gcc gct acco1 acsi1 ac_t at gtt ; run; quit;   Note that in the nonlinear AIDS model, two additional parameters are involved, namely, gtt and at. These parameters appear in the nonlinear price index specification. The parameter estimates for the nonlinear AIDS model are given in Figure 6 .   The MODEL Procedure Nonlinear ITSUR Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| Label gbb 0.037238 0.0176 2.11 0.0373   gbp -0.00538 0.00834 -0.65 0.5203   gbc -0.03931 0.0232 -1.69 0.0938   gbt 0.007449 0.0114 0.66 0.5137   gpp 0.043685 0.0111 3.94 0.0002   gpc -0.02431 0.0167 -1.46 0.1489   gpt -0.014 0.0112 -1.25 0.2140   gcc 0.071597 0.0195 3.67 0.0004   gct -0.00798 0.0191 -0.42 0.6766   gtt 0.014533 0.0215 0.68 0.5013   ab 0.377488 0.3819 0.99 0.3255   ap 0.158295 0.2667 0.59 0.5543   ac 0.632779 0.1670 3.79 0.0003   at -0.16856 0.2827 -0.60 0.5525   bb 0.02569 0.0423 0.61 0.5448   abco1 -0.01637 0.00162 -10.11 <.0001   absi1 -0.00224 0.00157 -1.43 0.1552   ab_t -0.00148 0.000049 -30.07 <.0001   bp 0.01255 0.0295 0.43 0.6716   apco1 0.00945 0.00117 8.10 <.0001   apsi1 0.006587 0.00112 5.86 <.0001   ap_t 0.000275 0.000036 7.69 <.0001   bc -0.05947 0.0185 -3.22 0.0018   acco1 -0.00554 0.000687 -8.06 <.0001   acsi1 -0.00094 0.000653 -1.44 0.1523   ac_t 0.000899 0.000026 34.15 <.0001   Restrict0 -170.721 74.3756 -2.30 0.0209 gbb + gbp + gbc + gbt = 0 Restrict1 13.41517 92.9286 0.14 0.8861 gbp + gpp + gpc + gpt = 0 Restrict2 282.4448 115.5 2.44 0.0137 gbc + gpc + gcc + gct = 0 Restrict3 -1.55944 0.7147 -2.18 0.0283 gbt + gpt + gct + gtt = 0 Restrict4 21.02431 6.7622 3.11 0.0015 ab + ap + ac + at = 1 Figure 6: Parameter Estimates for the Nonlinear AIDS Model Figure 7 illustrates the residuals obtained for the beef, pork, and chicken equations.     proc gplot data = resid2; plot w_beef*date w_pork*date w_chick*date / overlay cframe=ligr haxis=axis1 vaxis=axis2 vref=0; title 'Residual Plots'; footnote1 c=blue ' * beef ' c=red ' . pork ' c=green ' o chicken '; symbol1 c=blue i=join v=star; symbol2 c=red i=join v=dot; symbol3 c=green i = join v = circle ; axis1 label=('Time') ; axis2 label=(angle=90 'Residuals'); run; quit; Figure 7: Beef, Pork, and Chicken Residuals   The estimated parameters are stored in data sets and will be used in calculating price and income elasticities in the example "Calculating Elasticities in an Almost Ideal Demand System." Acknowledgment   The SAS program used in this example is based on code provided by Dr. Barry Goodwin.   References   Deaton, A., and Muellbauer, J. (1980a), "An Almost Ideal Demand System," American Economic Review, 70, 312-336. Deaton, A., and Muellbauer, J. (1980b), Economics and Consumer Behavior, Cambridge, England: Cambridge University Press. Piggott, N. (1997), The Benefits and Costs of Generic Advertising of Agricultural Commodities, Ph.D. dissertation, University of California, Davis. SAS Institute Inc. (1999), SAS/ETS User's Guide, Version 8, Cary, NC: SAS Institute Inc. USDA-ERS Electronic Data Archive, Red Meats Yearbook, housed at Cornell University's Mann Library, [http://usda.mannlib.cornell.edu/], accessed 25 September 2001. USDA-ERS Electronic Data Archive, Poultry Yearbook, housed at Cornell University's Mann Library, [http://usda.mannlib.cornell.edu/], accessed 25 September 2001. ... View more

Overview The process by which assets are priced is one of the fundamental problems in financial economics. One approach, developed by Lucus (1978), asserts that individuals hold assets in order to optimize their intertemporal consumption. In his classic model, individuals seek to maximize a time-additive intertemporal discounted utility function that depends upon stochastic consumption.   If an asset has a real rate of return Rt, an individual may be able to increase his utility by deferring consumption from the current period and investing in the asset in order to consume in a later period. The relative attractiveness between current and future consumption affects the asset's price as reflected in its return. An important implication is that changes in consumption should mirror changes in asset prices. Following Hansen and Singleton (1982), the agent maximizes the above utility function subject to the sequence of budget constraints where N is the number of assets in the economy, P it is the value of asset i in period t (price plus any cash payout), Q it is the amount of asset i owned by the individual at the end of period t, and W t is real labor income at date t. The first-order conditions yield the following Euler equations   where R t +1 is 1 plus the real rate of return on asset i, (P i,t+1 /P i,t ), and U' is marginal utility. This implies that, in equilibrium, the value of an asset is its discounted future payoff weighted by the trade-off between future and present consumption, conditioned on the public information set at time t, I t . By assuming a specific utility function you can relate the Euler equation that characterizes the agent's optimal solution to the returns that you see in the market. Assuming a power utility function [Note: this is a restrictive assumption. Refer to Pratt (1964) for a discussion on risk-aversion and utility functions.] with risk aversion coefficientα > 0, the above equation can be written in the estimable form where Z t is a vector of instruments that represent the public's information set at time t. This yields a set of moment conditions corresponding to the equilibrium asset prices and provides the proper environment to apply Hansen's (1982) generalized method of moments (GMM) to estimate the risk-aversion and time-preference parameters.   Analysis   With the assumption of a power utility function you can estimate a simple consumption-based asset pricing model from the Euler equations The left side of this expression can be thought of as an error term ut+1 that should have conditional mean 0, given the information set at time t, under the hypothesis of rational expectations. This defines a set of orthogonality conditions E( ut +1 | Z t ) = 0, implying E( ut +|Z t ) = 0. If Z t is any subset of the variables in the current information set, these orthogonality conditions can be exploited by using Hansen's (1982) generalized method of moments (GMM) to estimate the parameters of the model and to test their implications. Briefly, the GMM estimator is computed by minimizing the quadratic form q = m ' W -1 m where and W is the asymptotic variance/covariance matrix for the orthogonality conditions m. To perform the estimation requires data on some measure of consumption, returns for the assets of interest, and instruments from the public's information set. Of these, measures of consumption are the most difficult to obtain. Ferson and Harvey (1992) provide an informative discussion on the intricacies of consumption data. The data for this example are the same as those used by Ferson and Harvey (1992) in estimating equation (3) from their paper. They use quarterly observations from the second quarter of 1947 (1947.6) through the fourth quarter of 1987 (1987.12). Seasonally adjusted nondurables deflated by seasonally adjusted personal consumption deflators are used as the measure of consumption; the four assets used are the real returns from government and corporate bonds and the smallest and largest equity decile real returns. Real returns are nominal returns deflated by the price index corresponding to the consumption growth measure. In addition, four lags of the real Treasury bill return and four lags of real consumption growth, both deflated by the overall Con sumer Price Index (CPI), are used as the instruments, Z t . In the MODEL procedure, the endogenous variable is the percentage change in consumption, given by the variable CONRAT. The exogenous variables are the real returns for government bonds (GB), corporate bonds (CB), and the smallest and largest equity deciles (D1 and D10, respectively). One difference between the MODEL procedure and many other estimation tools is that the parameters to be estimated need to be specified explicitly. For this example, BETA denotes the time discount factor and ALPHA denotes the concavity parameter. Because computing parameter estimates for nonlinear equations requires an iterative process, you can supply the starting values. This example uses starting values of 1.0 for both parameters; the default starting values are 0.0001 for all parameters. proc model data=harvey; endogenous conrat; exogenous gb cb d1 d10 ; parms beta 1.0 alpha 1.0;   You can also transform variables for use in the estimation. The following statements create the four lagged values of real consumption growth (CINST) and the real Treasury bill return (RINST) to be used as instruments in the model. lc1 = lag(cinst); lc2 = lag2(cinst); lc3 = lag3(cinst); lc4 = lag4(cinst); ltb1 = lag(rinst); ltb2 = lag2(rinst); ltb3 = lag3(rinst); ltb4 = lag4(rinst);   The requisite moment conditions for the four assets are given by eq.h1 = beta * (1+conrat)**(-alpha) * (1+gb) - 1 ; eq.h2 = beta * (1+conrat)**(-alpha) * (1+cb) - 1 ; eq.h3 = beta * (1+conrat)**(-alpha) * (1+d1) - 1 ; eq.h4 = beta * (1+conrat)**(-alpha) * (1+d10) - 1 ;   You fit these equations by specifying the iterated GMM option using a Parzen kernel. Iterated GMM re-estimates the variance matrix at each iteration with the parameters determined by the GMM estimation from the previous iteration. See Ferson and Foerster (1994) for a discussion of iterated GMM. Iteration terminates when the variance matrix for the equation errors changes by less than the value of the CONVERGE= option. The default convergence value is 0.001. The KERNEL= option enables you to specify the smoothing function for the variance matrix. The INSTRUMENTS statement specifies the instrumental variables to be used in the ITGMM estimation method. This example uses the lagged variables created above. fit h1-h4 / itgmm kernel=(parzen,1,0); instruments lc1-lc4 ltb1-ltb4 ; run;     Consumption-Based Asset Pricing Model The MODEL Procedure Nonlinear ITGMM Summary of Residual Errors  Equation DF Model DF Error SSE MSE Root MSE R-Square Adj R-Sq h1 0.5 158.5 0.4277 0.00270 0.0519     h2 0.5 158.5 0.4111 0.00259 0.0509     h3 0.5 158.5 3.1048 0.0196 0.1400     h4 0.5 158.5 0.9861 0.00622 0.0789     Nonlinear ITGMM Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| beta 1.001594 0.00390 256.80 <.0001 alpha 0.645325 0.8546 0.76 0.4513 Number of Observations Statistics for System Used 159 Objective 0.2418 Missing 0 Objective*N 38.4454   For this simple model, the estimate of the time discount factor, BETA, is 1.00159 with an approximate standard error of 0.00390, while the estimate of the concavity parameter, ALPHA, is 0.64533 with an approximate standard error of 0.85462. In a time-separable model, such as this one, ALPHA can be thought of as the coefficient of relative risk aversion or as the inverse of the elasticity of consumption with respect to a fixed interest rate. The estimate of ALPHA in this example is not unusually large. These estimates are the same values reported by Ferson and Harvey (1992) in Panel B of Table IV in their paper, a four asset system with seasonally adjusted nondurables as the consumption measure. Any small differences can be attributed to the use of different kernel/bandwidth combinations. Let r be the number of unique instruments times the number of equations. The value r represents the number of orthogonality conditions imposed by the GMM method. Under the assumptions of the GMM method, r - p linearly independent combinations of the orthogonality conditions should be close to 0. The GMM estimates are computed by setting these combinations to 0. When r exceeds the number of parameters to be estimated, the OBJECTIVE*N, reported at the end of the estimation, is an asymptotically valid statistic to test the null hypothesis that the over-identifying restrictions of the model are valid. The OBJECTIVE*N is distributed as a chi-square with r - p degrees of freedom. The p-value of the test can be calculated with the following DATA step. data _null_; /* OBJECTIVE*N, degrees of freedom */ p = 1-probchi( 38.4454, 34 ); put p=; run;   In this example, OBJECTIVE*N = 38.445 is distributed as a chi-square with 34 degrees of freedom. This has an associated p-value of 0.275. Thus, you would not reject the null hypothesis that the overidentifying restrictions fit the model.   References   Ferson, W. and Foerster, S. (1994), ``Finite Sample Properties of the Generalized Method of Moments in Tests of Conditional Asset Pricing Models," The Journal of Financial Economics, 36, 29-55. Ferson, W. and Harvey, C. (1992), ``Seasonality and Consumption-Based Asset Pricing," The Journal of Finance, 47, 511-552. Hansen, L. (1982), ``Large Sample Properties of Generalized Method of Moments Estimators," Econometrica, 50, 1029-1084. Hansen, L. and Singleton, K. (1982), ``Generalized Instrumental Variables Estimation of Nonlinear Rational Expectations Models," Econometrica, 50, 1269-1286. Lucas, R. (1978), ``Asset Prices in an Exchange Economy," Econometrica, 48, 1149-1168. Pratt, J. (1964), ``Risk Aversion in the Small and in the Large," Econometrica, 32, 122-137. SAS Institute Inc. (1993), SAS/ETS User's Guide, Version 6, Second Edition, Cary, NC: SAS Institute Inc. ... View more

Overview   The generalized autoregressive conditional heteroscedasticity (GARCH) model of Bollerslev (1986) is an important type of time series model for heteroscedastic data. It explicitly models a time-varying conditional variance as a linear function of past squared residuals and of its past values. The GARCH process has been widely used to model economic and financial time-series data. Many extensions of the simple GARCH model have been developed in the literature. This example illustrates estimation of variants of GARCH models using the AUTOREG and MODEL procedures, which include the Simple GARCH model with normally distributed residuals GARCH model with t-distributed residuals GARCH model with Cauchy-distributed residuals GARCH model with generalized error distribution (GED) residuals GARCH-in-mean (GARCH-M) model Exponential GARCH (EGARCH) model Quadratic GARCH (QGARCH) model Glosten-Jagannathan-Runkle GARCH (GJR-GARCH) model Threshold GARCH (TGARCH) model Please note that parameter restrictions implied in the GARCH type models are not discussed in this example. If estimated parameters do not satisfy the desired restrictions in a specific model, the BOUNDS or RESTRICT statement can be used to explicitly impose the restrictions in PROC MODEL. For other examples of GARCH type models, see "Heteroscedastic Modeling of the Federal Funds Rate."   Details   The data used in this example are generated with the SAS DATA step. The following code generates a simple GARCH model with normally distributed residuals. %let df = 7.5; %let sig1 = 1; %let sig2 = 0.1 ; %let var2 = 2.5; %let nobs = 1000 ; %let nobs2 = 2000 ; %let arch0 = 0.1 ; %let arch1 = 0.2 ; %let garch1 = 0.75 ; %let intercept = 0.5 ; data normal; lu = &var2; lh = &var2; do i= -500 to &nobs ; /* GARCH(1,1) with normally distributed residuals */ h = &arch0 + &arch1*lu**2 + &garch1*lh; u = sqrt(h) * rannor(12345) ; y = &intercept + u; lu = u; lh = h; if i > 0 then output; end; run;     See the SAS program for more code that generates other types of GARCH models.   Simple GARCH Model with Normally Distributed Residuals   The simple GARCH(p,q) model can be expressed as follows. Let The residual   is modeled as   where   is i.i.d. with zero mean and unit variance, and where   is expressed as   In a standard GARCH model, is normally distributed. Alternative models can be specified by assuming different distributions for  , for example, the t distribution, Cauchy distribution, etc. To estimate a simple GARCH model, you can use the AUTOREG procedure. You use the GARCH= option to specify the GARCH model, and the (P= , Q= ) suboption to specify the orders of the GARCH model. proc autoreg data = normal ; /* Estimate GARCH(1,1) with normally distributed residuals with AUTOREG*/ model y = / garch = ( q=1,p=1 ) ; run ; quit ;   The AUTOREG procedure produces the following output given in Figure 1.1 for a GARCH model with normally distributed errors . Figure 1.1 Estimation Results using PROC AUTOREG The AUTOREG Procedure Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 0.4783 0.0559 8.55 <.0001 Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 0.4793 0.0323 14.86 <.0001 ARCH0 1 0.1159 0.0317 3.66 0.0003 ARCH1 1 0.2467 0.0389 6.35 <.0001 GARCH1 1 0.6972 0.0435 16.02 <.0001   You can also use the MODEL procedure to estimate a simple GARCH model. You must first specify the parameters in the model. Then specify the mean model and the variance model. The XLAG function returns the lag of the first argument if it is nonmissing. If the lag of the first argument is missing then the second argument is returned. The XLAG function makes it easy to specify the lag initialization for a GARCH process. The mse.y variable contains the value of the mean squared error for y at each iteration. These values are obtained automatically from first stage estimates, and are used to specify lagged values in estimation. /* Estimate GARCH(1,1) with normally distributed residuals with MODEL*/ proc model data = normal ; parms arch0 .1 arch1 .2 garch1 .75 ; /* mean model */ y = intercept ; /* variance model */ h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y) ; /* fit the model */ fit y / method = marquardt fiml ; run ; quit ;   Figure 1.2  shows the parameter estimates obtained by using the MODEL procedure. Figure 1.2 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.479341 0.0319 15.02 <.0001 arch0 0.115242 0.0345 3.34 0.0009 arch1 0.246811 0.0432 5.72 <.0001 garch1 0.697988 0.0494 14.13 <.0001   GARCH Model with t-Distributed Residuals To estimate a GARCH model with -distributed errors, you can use the AUTOREG procedure. You specify the GARCH(p,q) process with the GARCH=(p=,q=) option, and specify the distributed error structure with the DIST= option. /* Estimate GARCH(1,1) with t-distributed residuals with AUTOREG*/ proc autoreg data = t ; model y = / garch=( q=1, p=1 ) dist = t ; run ; quit;   The AUTOREG procedure produces the following output given in Figure 1.3  for a GARCH model with-distributed residuals. Figure 1.3 Estimation Results using PROC AUTOREG The AUTOREG Procedure Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 0.3997 0.0998 4.00 <.0001 Variable DF Estimate Standard Error t Value Approx Pr > |t| Variable Label Intercept 1 0.5294 0.0374 14.14 <.0001   ARCH0 1 0.1244 0.0328 3.80 0.0001   ARCH1 1 0.2919 0.0309 9.44 <.0001   GARCH1 1 0.7419 0.0199 37.27 <.0001   TDFI 1 0.1351 0.0250 5.41 <.0001 Inverse of t DF   You can also use the MODEL procedure to estimate the GARCH model with -distributed residuals. Use the ERRORMODEL statement to specify the -distributed residuals. First specify the dependent variable name, a tilde(~), and then the name of the error distribution with its parameters. The degrees of freedom for the distribution are also estimated as a parameter in the MODEL procedure. /* Estimate GARCH(1,1) with t-distributed residuals with MODEL*/ proc model data = t ; parms df 7.5 arch0 .1 arch1 .2 garch1 .75 ; /* mean model */ y = intercept ; /* variance model */ h.y = arch0 + arch1 * xlag(resid.y **2, mse.y) + garch1*xlag(h.y, mse.y); /* specify error distribution */ errormodel y ~ t(h.y,df); /* fit the model */ fit y / method=marquardt; run; quit; The MODEL procedure produces the following output given in Figure 1.4 for the GARCH model with -distributed errors. Figure 1.4 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear Liklhood Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.528964 0.0376 14.05 <.0001 df 7.468812 1.2402 6.02 <.0001 arch0 0.091711 0.0249 3.69 0.0002 arch1 0.215409 0.0248 8.69 <.0001 garch1 0.740497 0.0222 33.34 <.0001   Note that the MODEL procedure outputs the estimate of the degree of freedom for the distribution, while the AUTOREG procedure outputs the reciprocal of the estimated degree of freedom.   GARCH Model with Cauchy Distributed Residuals   You can also estimate a GARCH model with Cauchy-distributed errors in the MODEL procedure. The log likelihood function for GARCH with Cauchy-distributed residuals can be expressed as   We can then write the code using the model procedure. /* Estimate GARCH(1,1) with Cauchy distributed residuals */ proc model data = cauchy ; parms arch0 .1 arch1 .2 garch1 .75 intercept .5 ; mse_y = &var2 ; /* mean model */ y = intercept ; /* variance model */ h.y = arch0 + arch1 * xlag(resid.y ** 2, mse_y) + garch1 * xlag(h.y, mse_y); /* specify error distribution */ obj = log(h.y/((h.y**2+resid.y**2) * constant('pi'))); obj = -obj ; errormodel y ~ general(obj); /* fit the model */ fit y / method=marquardt; run; quit;   The MODEL procedure produces the following output given in Figure 1.5 for a GARCH fit with Cauchy distributed residuals. Figure 1.5 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear Liklhood Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.499152 0.00385 129.73 <.0001 arch0 0.024821 0.00350 7.10 <.0001 arch1 0.005344 0.00180 2.97 0.0030 garch1 0.678644 0.0417 16.29 <.0001   GARCH Model with Generalized Error Distribution Residuals You can also estimate a GARCH model with GED (generalized error distribution) residuals with the MODEL procedure. The log likelihood function for GARCH with GED residuals is expressed as where T is the sample size, Γ(•)is the gamma function, λ is a constant given by and   is a positive parameter governing the thickness of the tails of the distribution. Note that for   = 2, constant λ = 1, and the GED is the standard normal distribution. For more details about the generalized error distribution, see Hamilton (1994). To estimate a GARCH model with GED residuals, you can use the following code. /* Estimate GARCH(1,1) with generalized error distribution residuals */ proc model data = normal ; parms nu 2 arch0 .1 arch1 .2 garch1 .75; control mse.y = &var2 ; /*defined in data generating step*/ /* mean model */ y = intercept ; /* variance model */ h.y = arch0 + arch1 * xlag(resid.y ** 2, mse.y) + garch1 * xlag(h.y, mse.y); /* specify error distribution */ lambda = sqrt(2**(-2/nu)*gamma(1/nu)/gamma(3/nu)) ; obj = log(nu/lambda) -(1 + 1/nu)*log(2) - lgamma(1/nu)- .5*abs(resid.y/lambda/sqrt(h.y))**nu - .5*log(h.y) ; obj = -obj ; errormodel y ~ general(obj,nu); /* fit the model */ fit y / method=marquardt; run; quit;   The MODEL procedure produces the following output given in Figure 1.6  for a GARCH fit with residuals following a generalized error distribution. Figure 1.6 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear Liklhood Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.47119 0.0321 14.69 <.0001 nu 1.964288 0.1336 14.70 <.0001 arch0 0.172375 0.0359 4.80 <.0001 arch1 0.276887 0.0441 6.28 <.0001 garch1 0.637912 0.0476 13.39 <.0001   GARCH-M Model   Another type of GARCH model is the GARCH-M model, which adds the heteroscedasticity term directly into the mean equation. In this example, consider the following specification:   The residual   is modeled as where   is i.i.d. with zero mean and unit variance, and where   is expressed as The AUTOREG procedure enables you to specify the GARCH-M model with the MEAN= sub-option of the GARCH= option. The MEAN= option specifies the functional form of the GARCH-M model. The values of the MEAN= option are LINEAR, specifies the linear function LOG, specifies the log function SQRT, specifies the square-root function In this example, the square-root specification is considered. /* Estimate GARCH-M Model with PROC AUTOREG */ proc autoreg data= garchm ; model y = / garch=( p=1, q=1, mean = sqrt); run; quit;   The AUTOREG procedure produces the following output in Figure 1.7 for a GARCH-M model. Figure 1.7 Estimation Results using PROC AUTOREG The AUTOREG Procedure Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 1.2190 0.0480 25.39 <.0001 Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 0.4987 0.1431 3.49 0.0005 ARCH0 1 0.0831 0.0275 3.02 0.0025 ARCH1 1 0.1695 0.0289 5.86 <.0001 GARCH1 1 0.7914 0.0321 24.62 <.0001 DELTA 1 0.5206 0.1196 4.35 <.0001   You can also use the MODEL procedure to estimate the GARCH-M model. /* Estimate GARCH-M Model with PROC MODEL */ proc model data = garchm ; parms arch0 .1 arch1 .2 garch1 .75 gamma .5 ; h = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y); y = intercept + gamma*sqrt(h) ; h.y = h ; fit y / fiml method = marquardt; run; quit;   This PROC MODEL step produces the following output as shown in Figure 1.8 . Figure 1.8 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| arch0 0.082298 0.0268 3.07 0.0022 arch1 0.17043 0.0284 6.00 <.0001 garch1 0.79191 0.0318 24.87 <.0001 gamma 0.516689 0.1225 4.22 <.0001 intercept 0.502515 0.1488 3.38 0.0008     EGARCH Model   The Exponential GARCH (EGARCH) model was proposed by Nelson (1991). It models the conditional variance of   as follows: where The AUTOREG procedure also supports the EGARCH model. You can use the TYPE=EXP suboption of the GARCH= option to specify the EGARCH model. /* Estimate EGARCH Model with PROC AUTOREG */ proc autoreg data= egarch ; model y = / garch=( q=1, p=1 , type = exp) ; run; quit;   This produces the following output as shown in Figure 1.9 . Figure 1.9 Estimation Results using PROC AUTOREG The AUTOREG Procedure Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 0.4790 0.0381 12.59 <.0001 Variable DF Estimate Standard Error t Value Approx Pr > |t| Intercept 1 0.4854 0.0361 13.43 <.0001 EARCH0 1 0.0960 0.0366 2.63 0.0087 EARCH1 1 0.2322 0.0744 3.12 0.0018 EGARCH1 1 0.7206 0.0961 7.50 <.0001 THETA 1 0.4403 0.1806 2.44 0.0148   You can also estimate the EGARCH model using the MODEL procedure.   /* Estimate EGARCH Model with PROC MODEL */ proc model data = egarch ; parms earch0 .1 earch1 .2 egarch1 .75 theta .65 ; /* mean model */ y = intercept ; /* variance model */ if (_obs_ = 1 ) then h.y = exp( earch0 + egarch1 * log(mse.y) ); else h.y = exp(earch0 + earch1*zlag(g) + egarch1*log(zlag(h.y))) ; g = theta*(-nresid.y) + abs(-nresid.y) - sqrt(2/constant('pi')) ; /* fit the model */ fit y / fiml method = marquardt ; run; quit;   Note that in this example, the parameter   is set to be equal to 1. The nresid.y variable gives the normalized residual of the variable y, i.e.,  . Note also that   if  . The PROC MODEL code produces the following output for the EGARCH model as shown in Figure 1.10 . Figure 1.10 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.485374 0.0368 13.18 <.0001 earch0 0.095713 0.0361 2.65 0.0081 earch1 0.233174 0.0690 3.38 0.0007 egarch1 0.721372 0.0934 7.73 <.0001 theta 0.434851 0.1790 2.43 0.0153     QGARCH Model   The Quadratic GARCH model was proposed by Sentana (1995) to model asymmetric effects of positive and negative shocks. A simple Quadratic GARCH(1,1) model describes the residual process   as where  is i.i.d. with zero mean and unit variance, and   where φ is the asymmetric parameter that helps to separately identify the impact of positive and negative shocks on volatility. The following code estimates a simple Quadratic GARCH(1,1) model in the MODEL procedure. /* Estimate Quadratic GARCH (QGARCH) Model */ proc model data = qgarch ; parms arch0 .1 arch1 .2 garch1 .75 phi .2; /* mean model */ y = intercept ; /* variance model */ h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y) + phi*xlag(-resid.y,mse.y); /* fit the model */ fit y / method = marquardt fiml ; run ; quit ;   Note that in specifying the equation for , you need to add a negative sign in front of the residual term, since PROC MODEL gives the negative of the residuals. This also applies to the GJR-GARCH model and the TGARCH model, to be discussed later in the example. The code produces the following output shown in Figure 1.11 for a Quadratic GARCH(1,1) model. Figure 1.11 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear FIML Summary of Residual Errors Equation DF Model DF Error SSE MSE Root MSE R-Square Adj R-Sq y 5 995 2090.0 2.1005 1.4493 -0.0004 -0.0045 resid.y   995 996.5 1.0015 1.0007     Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.510972 0.0356 14.37 <.0001 arch0 0.108659 0.0281 3.86 0.0001 arch1 0.166853 0.0294 5.67 <.0001 garch1 0.774571 0.0341 22.69 <.0001 phi 0.16391 0.0426 3.84 0.0001     GJR-GARCH Model   Another asymmetric GARCH process is the GJR-GARCH model of Glosten, Jagannathan and Runkle (1993). They propose modeling  , where   is i.i.d. with zero mean and unit variance, and where   if   and   if  . You can use the following code to estimate a GJR-GARCH(1,1) model. /* Estimate GJR-GARCH Model */ proc model data = gjrgarch ; parms arch0 .1 arch1 .2 garch1 .75 phi .1; /* mean model */ y = intercept ; /* variance model */ if zlag(resid.y) > 0 then h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y) ; else h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y) + phi*xlag(resid.y**2,mse.y) ; /* fit the model */ fit y / method = marquardt fiml ; run ; quit ;   This produces the following output shown in Figure 1.12  with the GJR-GARCH model. Figure 1.12 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear FIML Summary of Residual Errors Equation DF Model DF Error SSE MSE Root MSE R-Square Adj R-Sq y 5 995 13387.5 13.4548 3.6681 -0.0000 -0.0040 resid.y   995 996.2 1.0012 1.0006     Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.557155 0.0513 10.85 <.0001 arch0 0.099621 0.0348 2.87 0.0042 arch1 0.185367 0.0386 4.80 <.0001 garch1 0.771999 0.0269 28.71 <.0001 phi 0.089578 0.0509 1.76 0.0788     TGARCH Model   The Threshold GARCH model (TGARCH) of Zakoian (1994) is similar to the GJR GARCH, but it specifies the conditional standard deviation instead of conditional variance: where   if  , and   if  . Similarly,   if  , and   if  . You can estimate TGARCH(1,1) model using the following code. /* Estimate Threshold Garch (TGARCH) Model */ proc model data = tgarch ; parms arch0 .1 arch1_plus .1 arch1_minus .1 garch1 .75 ; /* mean model */ y = intercept ; /* variance model */ if zlag(resid.y) < 0 then h.y = (arch0 + arch1_plus*zlag(-resid.y) + garch1*zlag(sqrt(h.y)))**2 ; else h.y = (arch0 + arch1_minus*zlag(-resid.y) + garch1*zlag(sqrt(h.y)))**2 ; /* fit the model */ fit y / method = marquardt fiml ; run ; quit ;   This produces the following output for the TGARCH model as shown in Figure 1.13. Figure 1.13 Estimation Results using PROC MODEL The MODEL Procedure Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| intercept 0.504727 0.0124 40.57 <.0001 arch0 0.160549 0.0456 3.52 0.0004 arch1_plus 0.076379 0.0361 2.11 0.0347 arch1_minus 0.14484 0.0307 4.72 <.0001 garch1 0.632826 0.1161 5.45 <.0001     References   Bollerslev, T. (1986), "Generalized Autoregressive Conditional Heteroskedasticity," Journal of Econometrics, 31, 307-327. Glosten, L., Jagannathan, R. and Runkle, D. (1993), "On the Relation between the Expected Value and the Volatility of the Nominal Excess Return on Stocks," Journal of Finance, 48(5), 1779-1801. Hamilton, J. D. (1994), Time Series Analysis, Princeton, NJ: Princeton University Press. Nelson, B. (1991), "Conditional Heteroskedasticity in Asset Returns: A New Approach," Econometrica, 59, 347-370. SAS Institute Inc. (1999), SAS/ETS User’s Guide, Version 8, Cary, NC: SAS Institute Inc. Sentana, E. (1995), "Quadratic ARCH Models," Review of Economic Studies, 62, 639-661. Zakoian, M. (1994), "Threshold Heteroscedastic Models," Journal of Economic Dynamics and Control, 18, 931-955. ... View more

Overview   The Capital Asset Pricing  Model (CAPM) is one of the most common methods of relating the sensitivity of an individual company's stock return to the return of the market as a whole. Specifically, the CAPM linear relationship can be written as where r j , r f , and r m are returns to security j, the risk-free asset and the market return, respectively. The term is the ratio of the standard deviations of returns on security j to the market. A full development of the CAPM can be found in any introductory finance textbook (for example, Brealey and Myers 1988), and so is omitted here. There are, however, several econometric issues involved in the more general model where   is an intercept term, and   is an IID stochastic disturbance term with mean zero and variance  . Two straightforward tests of this model are a test of the null against the alternative that it is not equal to 0. A rejection would cast some doubt on the theoretical model as α j is assumed to be 0 in the first equation. a test of the null  against the alternative that it is not equal to 1. Failure to reject would imply that the movement of asset prices of a particular company is the same as that of the market as a whole. This example uses stock returns from the Tandy Corporation to estimate a CAPM and perform some general tests of the model. The CAPM relates the sensitivity of an individual company's stock returns to the returns of the market as a whole. Estimating a model for a particular firm requires data on the market rate of return (typically a composite index such as the S&P 500 or the Dow Jones Industrial Average), the risk-free rate of return (usually a short-term Treasury bill), and stock returns from the company of interest. The data for this example consist of monthly observations from January 1978 through December 1987 on the market return, the risk-free rate, and the return on the Tandy Corporation's common stock. Using a SAS DATA step, enter the data and create two new variables, R_TANDY = TANDY - R_F and R_MKT = R_M - R_F, that correspond to the risk premiums for the Tandy Corporation and the Market. A risk premium is the excess return of a security over the risk-free rate or, rather, the extra return that investors require for bearing risk. data tandy; input r_m r_f tandy @@; r_tandy = tandy - r_f; r_mkt = r_m -r_f; label r_m='Market Rate of Return' r_f='Risk-Free Rate of Return' tandy='Rate of Return for Tandy Corporation' r_tandy='Risk Premium for Tandy Corporation' r_mkt='Risk Premium for Market'; datalines; -.045 .00487 -.075 .063 .00491 .055 .071 .00528 .194 ... more data ... ; run;   It is always a good idea to look at the data to be sure that it was input as expected. A simple plot allows for a visual assessment of the appropriateness of the proposed model. The following code uses the GPLOT procedure to plot the risk premium of the Tandy Corporation against the risk premium of the Market. proc gplot data=tandy; plot r_tandy * r_mkt / haxis=axis1 hminor=4 cframe=ligr vaxis=axis2 vminor=4; symbol1 c=blue v=star; axis1 order=(-0.3 to 0.2 by 0.1); axis2 label=(angle=90 'Tandy Corp. Risk Premium') order=(-0.4 to 0.6 by 0.2); title 'TANDY Corporation CAPM Example'; title2'Plot of Risk Premiums'; title3'Tandy Corporation versus the Market'; run;     The plot of the Tandy data set above indicates that there is a positive relationship between the risk premium of the market and the risk premium of the Tandy Corporation's stock. This positive relationship seems sensible: as the market return increases, so should the return of most stocks. Based on this visual evidence, a linear model between R_TANDY and R_MKT is not unreasonable. The AUTOREG procedure specifies a linear regression of R_TANDY on R_MKT. Note that a constant term is automatically included in the model unless the NOINT option is specified. The DWPROB option requests the Durbin-Watson test for autocorrelation of the residuals. The TEST statement enables you to perform hypothesis tests on parameter estimates; the following TEST statement is requesting a test that the coefficient of R_MKT is equal to 1. In addition, the OUTPUT data set TANDYOUT contains the predicted values of R_TANDY, the residuals of the regression, and the upper and lower bounds of the 90% confidence interval. proc autoreg data=tandy; model r_tandy = r_mkt / dwprob; test r_mkt = 1; output out=tandyout p=p r=r ucl=u lcl=l alphacli=.10; title2; title3; run;   The AUTOREG Procedure Dependent Variable r_tandy   Risk Premium for Tandy Corporation    Ordinary Least Squares Estimates SSE 1.31740165 DFE 118 MSE 0.01116 Root MSE 0.10566 SBC -191.29942 AIC -196.8744 Regress R-Square 0.3191 Total R-Square 0.3191 Durbin-Watson 1.9334 Pr < DW 0.3587 Pr > DW 0.6413     Variable DF Estimate Standard Error t Value Approx Pr > |t| Variable Label Intercept 1 0.0107 0.009698 1.10 0.2740   r_mkt 1 1.0500 0.1412 7.44 <.0001 Risk Premium for Market Test 1 Source DF Mean Square F Value Pr > F Numerator 1 0.001400 0.13 0.7239 Denominator 118 0.011164     The output provides information about the appropriateness of the model chosen. The parameter estimate for , the intercept, is not significantly different from 0. The estimate for β, the coefficient of R_MKT, is 1.05 and is highly significant. So, not altogether unexpectedly, it seems that the market risk premium is a factor in the risk premium for the Tandy Corporation. The F-statistic for the test that the coefficient of R_MKT is equal to one is 0.1254 with an associated p-value of 0.724, implying that there is little evidence that this coefficient is different from 1.0. Thus there is little evidence that the returns of the Tandy stock were any more or less volatile than the returns of the market. R2 measures the proportion of the total variation of Y explained by the regression of Y on the independent variable, X. The R2 value of 0.319 means that about 32% of the variation in the risk premium of Tandy Corporation stock can be explained by the risk premium of the market. Or, put another way, 68% of the risk premium of Tandy stock is firm specific. The Durbin-Watson d statistic tests for autocorrelation and lack of independence of residuals, which is a common problem in time series data. The d statistic ranges from 0 to 4. A value close to 2.0 indicates that you cannot reject the null hypothesis of no autocorrelation. The calculated value of the d statistic in the CAPM for the Tandy Corporation is 1.933, and it has a p-value of 0.359. This value indicates that, for this model, you cannot reject the null hypothesis of independent (nonautocorrelated) residuals at typically chosen levels of significance. Using PROC GPLOT, you can investigate the residuals for patterns that would indicate a violation of the underlying assumptions of the model. proc gplot data=tandyout; plot r * r_mkt / haxis=axis1 hminor=4 cframe=ligr vaxis=axis2 vminor=4 vref=0.0; symbol1 c=green v=star; axis1 order=(-0.3 to 0.2 by 0.1); axis2 label=(angle=90 'Tandy Corp. Risk Premium') order=(-0.3 to 0.4 by 0.1); title 'TANDY Corporation CAPM Example'; title2'OLS Residuals versus Market Risk Premium'; run;   This residual plot does not reveal any obvious patterns; the residuals appear randomly distributed about the reference line.   References Berndt, Ernst R. (1991), The Practice of Econometrics: Classic and Contemporary, New York: Addison-Wesley. Brealey, Richard A. and Myers, Stewart C. (1988), Principles of Corporate Finance, Third Edition, New York: McGraw-Hill. SAS Institute Inc. (1993), SAS/ETS Software: Applications Guide 2, Version 6, First Edition: Econometric Modeling, Simulation, and Forecasting, Cary, NC: SAS Institute Inc., 15-30. ... View more

Overview   Many economic and business variables are affected by seasonal factors. For example, power usage is highest in the months when temperatures are most extreme. The most common type of seasonality is variation due to the time of year, but other types of seasonality are also found in time series data. Seasonal models are often multiplicative rather than additive. A multiplicative model includes the product of one or more nonseasonal parameters with one or more seasonal parameters. For example, a multiplicative model with both autoregressive and moving average terms (an ARMA model) and with yearly seasonality for a time series, y t , can be written as: where is the intercept parameter.  is the nonseasonal first-order autoregressive parameter. is the seasonal autoregressive parameter.  is the nonseasonal first-order moving average parameter. is the seasonal moving average parameter. To identify a seasonal model, you need to examine the autocorrelation function (ACF) and the inverse autocorrelation function (IACF) plots. For multiplicative MA processes, there are small spikes in the ACF plot q lags before and after the seasonal lag, where q is the number of nonseasonal MA parameters necessary to model the data. These small spikes are usually in the opposite direction of the seasonal spike. For example, a multiplicative MA(1, 12) process typically has small spikes at lags 11 and 13 on either side of, and in the opposite direction of, a large spike at lag 12. An additive MA process typically has small spikes q lags before the seasonal lag, where q is the number of nonseasonal MA parameters necessary to model the data. For example, an additive MA(1, 12) process typically has a small spike at lag 11 and a larger spike at lag 12. To identify an AR process, look for the patterns described previously in the IACF plot rather than in the ACF plot. If a process contains both AR and MA components, the patterns may appear in both the ACF and IACF plots. This example develops an ARMA model for steel shipments from U.S. steel mills.   Analysis The identification and estimation of Autoregressive Integrated Moving Average (ARIMA) models is more of an art than a science. Generally, the most parsimonious model fitting the data is considered the best. This example uses steel shipments data taken from Metal Statistics 1993. The values represent monthly totals of steel products shipped from U.S. steel mills, in thousands of net tons, for the period from January 1984 to December 1991. The following statements create the data set STEEL.   data steel; input date:monyy5. steelshp @@; format date monyy5.; title 'U.S. Steel Shipments Data'; title2 '(thousands of net tons)'; datalines; JAN84 5980 FEB84 6150 MAR84 7240 APR84 6472 MAY84 6948 JUN84 6686 JUL84 5820 AUG84 6033 SEP84 5454 OCT84 6087 NOV84 5317 DEC84 4867 ... more data lines ... ;   The analysis performed by the ARIMA procedure is divided into three stages, corresponding to the stages described by Box and Jenkins (1976). The IDENTIFY, ESTIMATE, and FORECAST statements perform these three stages. In the identification stage, you use the IDENTIFY statement to specify the response series and identify candidate ARIMA models for it. The IDENTIFY statement reads time series that are to be used in later statements, possibly differencing them, and computes autocorrelations, inverse autocorrelations, partial autocorrelations, and cross correlations. The analysis of this output usually suggests one or more ARIMA models that could be fit. The VAR= option specifies the variable to be identified. proc arima data=steel; i var=steelshp; run;   U.S. Steel Shipments Data (thousands of net tons) The ARIMA Procedure Autocorrelations Lag Covariance Correlation -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1 Std Error 0 406442 1.00000 |                    |********************| 0 1 262630 0.64617 |                .   |*************       | 0.102062 2 261597 0.64363 |              .     |*************       | 0.138258 3 235909 0.58042 |             .      |************        | 0.166570 4 168515 0.41461 |             .      |********            | 0.186451 5 201896 0.49674 |            .       |**********          | 0.195820 6 129000 0.31739 |            .       |****** .            | 0.208533 7 152701 0.37570 |           .        |********.           | 0.213506 8 113117 0.27831 |           .        |******  .           | 0.220285 9 127532 0.31378 |           .        |******  .           | 0.223918 10 137000 0.33707 |           .        |******* .           | 0.228452 11 130723 0.32163 |           .        |******  .           | 0.233575 12 200408 0.49308 |          .         |**********          | 0.238144 13 112496 0.27678 |          .         |******   .          | 0.248551 14 135119 0.33244 |          .         |*******  .          | 0.251741 15 103295 0.25414 |          .         |*****    .          | 0.256273 16 62982.090 0.15496 |          .         |***      .          | 0.258885 17 108381 0.26666 |          .         |*****    .          | 0.259850 18 42836.479 0.10539 |         .          |**        .         | 0.262685 19 65840.039 0.16199 |         .          |***       .         | 0.263125 20 37765.859 0.09292 |         .          |**        .         | 0.264162 21 27790.106 0.06837 |         .          |*         .         | 0.264502 22 40303.846 0.09916 |         .          |**        .         | 0.264686 23 46097.710 0.11342 |         .          |**        .         | 0.265073 24 76317.464 0.18777 |         .          |****      .         | 0.265578 "." marks two standard errors The large spike at lag 12 in the ACF plot provides evidence that the steel shipments time series has a seasonal autoregressive component. The lack of a large spike at lag 24 indicates that the series is stationary at the seasonal level. U.S. Steel Shipments Data (thousands of net tons) The ARIMA Procedure Inverse Autocorrelations Lag Correlation -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1 1 -0.37291 |             *******|   .                | 2 0.08136 |                .   |** .                | 3 -0.31032 |              ******|   .                | 4 0.16197 |                .   |***.                | 5 -0.20750 |                ****|   .                | 6 0.16115 |                .   |***.                | 7 -0.02341 |                .   |   .                | 8 0.06910 |                .   |*  .                | 9 0.00628 |                .   |   .                | 10 0.02046 |                .   |   .                | 11 0.02875 |                .   |*  .                | 12 -0.23279 |               *****|   .                | 13 0.03755 |                .   |*  .                | 14 0.04050 |                .   |*  .                | 15 0.03498 |                .   |*  .                | 16 0.09969 |                .   |** .                | 17 -0.10703 |                . **|   .                | 18 0.04901 |                .   |*  .                | 19 -0.08634 |                . **|   .                | 20 0.02281 |                .   |   .                | 21 0.00844 |                .   |   .                | 22 0.10510 |                .   |** .                | 23 -0.10923 |                . **|   .                | 24 0.02676 |                .   |*  .                | The spikes at lags 1 and 3 in the IACF plot indicate that other components are necessary to fit an adequate model. The null hypothesis of white noise residuals is resoundingly rejected. U.S. Steel Shipments Data (thousands of net tons) The ARIMA Procedure Autocorrelation Check for White Noise To Lag Chi-Square DF Pr > ChiSq Autocorrelations 6 170.51 6 <.0001 0.646 0.644 0.580 0.415 0.497 0.317 12 255.47 12 <.0001 0.376 0.278 0.314 0.337 0.322 0.493 18 296.96 18 <.0001 0.277 0.332 0.254 0.155 0.267 0.105 24 309.34 24 <.0001 0.162 0.093 0.068 0.099 0.113 0.188 In the estimation and diagnostic checking stage, you use the ESTIMATE statement to specify the ARIMA model to fit to the variable specified in the previous IDENTIFY statement and to estimate the parameters of that model. The ESTIMATE statement also produces diagnostic statistics to help you judge the adequacy of the model. Significance tests for parameter estimates indicate whether some terms in the model may be unnecessary. Goodness-of-fit statistics aid in comparing this model to others. Tests for white noise residuals indicate whether the residual series contains additional information that might be used by a more complex model. If the diagnostic tests indicate problems with the model, you try another model, then repeat the estimation and diagnostic checking stage. The following statement fits a seasonal ARMA model to the time series. In the syntax of the ESTIMATE statement, the two multiplicative AR terms, denoted by the P= option, are enclosed in separate parentheses. The two additive MA terms, denoted by the Q= option, are separated by a space within a single set of parentheses. e p=(2)(12) q=(1 3); run;   U.S. Steel Shipments Data (thousands of net tons) The ARIMA Procedure Autocorrelation Check of Residuals To Lag Chi-Square DF Pr > ChiSq Autocorrelations 6 2.42 2 0.2979 -0.009 -0.051 0.071 0.070 0.104 0.018 12 3.63 8 0.8891 -0.084 0.032 -0.024 0.013 -0.033 -0.035 18 11.86 14 0.6176 -0.082 0.168 0.014 -0.137 0.107 0.073 24 16.16 20 0.7066 0.023 0.019 -0.010 -0.047 0.174 -0.000 Model for variable steelshp Estimated Mean 6057.122 Autoregressive Factors Factor 1: 1 - 0.54234 B**(2) Factor 2: 1 - 0.64802 B**(12) Moving Average Factors Factor 1: 1 + 0.55505 B**(1) + 0.43689 B**(3) The Autocorrelation Check of Residuals shows that none of the Q-statistics are statistically significant. This indicates that the model provides an adequate fit to the data. U.S. Steel Shipments Data (thousands of net tons) The ARIMA Procedure Conditional Least Squares Estimation Parameter Estimate Standard Error t Value Approx Pr > |t| Lag MU 6057.1 232.96713 26.00 <.0001 0 MA1,1 -0.55505 0.08021 -6.92 <.0001 1 MA1,2 -0.43689 0.07936 -5.51 <.0001 3 AR1,1 0.54234 0.09903 5.48 <.0001 2 AR2,1 0.64802 0.09392 6.90 <.0001 12 Constant Estimate 975.7391 Variance Estimate 126334.1 Std Error Estimate 355.4351 AIC 1404.983 SBC 1417.805 Number of Residuals 96 All of the estimated parameters have relatively large t-statistics, which indicates that these parameters cannot be omitted from the model. In the forecasting stage, you use the FORECAST statement to forecast future values of the time series and to generate confidence intervals for these forecasts from the ARIMA model produced by the preceding ESTIMATE statement. The following statements produce forecasts and upper and lower 95% confidence limits for 12 future periods and creates the output data set STEEL2. f lead=12 out=steel2 id=date interval=month noprint; run;   To prepare the output data set for plotting, change the values for the forecasts and confidence limits to missing for all dates prior to the future forecast periods. data steel3; set steel2; if date lt '01jan92'd then do; forecast=.; l95=.; u95=.; end; run;   Use the GPLOT procedure to plot the data. proc gplot data=steel3; format date year4.; plot steelshp*date=1 forecast*date=2 l95*date=3 u95*date=3 / overlay cframe=ligr haxis=axis1 vaxis=axis2 vminor=1 href='01jan92'd; title 'U.S. Steel Shipments Data'; title2 '(thousands of net tons)'; axis1 offset=(1 cm) label=('Year') minor=none order=('01jan84'd to '01jan93'd by year); axis2 label=(angle=90 'Steel Shipments') order=(4500 to 8500 by 1000); symbol1 c=blue i=join l=1 v=star; symbol2 c=red i=join l=1 v=F; symbol3 c=green i=join l=20; run; quit;     The values of the original steel shipments time series are plotted with the star symbol. The forecasts are plotted with the F symbol, and the upper and lower 95% confidence limits for the forecasts are plotted with dashed lines. Because the model fit to the steel shipments data includes a seasonal component, the forecasts do not follow a simple linear trend. Instead, the forecasts show variability due to the season (month of the year).   References Box, G.E.P. and Jenkins, G.M. (1976), Time Series Analysis: Forecasting and Control, San Francisco: Holden-Day. Chilton Publications (1993), Metal Statistics 1993, New York: Chilton Publications. Hamilton, J. (1994), Time Series Analysis, Princeton, NJ: Princeton University Press. SAS Institute Inc. (1996), Forecasting Examples for Business and Economics Using the SAS System, Cary, NC: SAS Institute Inc. SAS Institute Inc. (1993), SAS/ETS User's Guide, Version 6, Second Edition, Cary, NC: SAS Institute Inc. ... View more

Overview Figure 1.1: Heteroscedastic residuals Modeling non-constant variance, or heteroscedasticity, improves the efficiency of estimates of the parameters associated with the mean of a series and provides insight into the volatility of a series. SAS/ETS software provides capability to do linear and nonlinear regression with heteroscedastic models using the AUTOREG and MODEL procedures, respectively. This example makes use of the MODEL procedure.   Details One of the key assumptions of regression analysis is that the variance of the errors is constant across observations. This assumption is often violated when modeling time series or panel data, resulting in inefficient parameter estimates and inaccurate forecast error variance. Consider the following general model: where   For models that are homoscedastic, h(t)=1. If you have a model which is heteroscedastic with known form, you can improve the efficiency of the estimates by performing a weighted regression. The weight variable, using this notation, is If the errors for a model are heteroscedastic and the functional form of the variance is known, the model for the variance can be estimated along with the regression function. To specify a functional form for the variance, assign the function to an H.var variable where var is the equation variable. For example, if you want to estimate the scale parameter for the variance of a simple OLS regression model   y = a + b * x   you can specify   proc model data=s; y = a + b * x; h.y = sigma**2; fit y;   The above model has a constant variance, sigma**2. PROC MODEL estimates the constant-variance model by default, therefore the h.y equation is unnecessary in this case. Now, consider the same model with the following functional form for the variance: This model is written as proc model data=s; y = a + b * x; h.y = sigma**2 * x**alpha; fit y;   In addition to estimating the parameters, a and b, the MODEL procedure also estimates the parameters sigma and alpha, which are associated with the error variance. The Federal Funds Rate     Figure 1.2: The Weekly Federal Funds Rate The following example uses weekly effective Federal funds data created from an average of daily figures from January 1961 to March 2000. A simple, commonly used rate model is the Vasicek model: The following PROC MODEL statements are used to fit this model: proc model data=weeklyrates; ffrate = lag(ffrate) + kappa * (theta - lag(ffrate)); lag_ffrate = lag( ffrate ); label kappa = "Speed of Mean Reversion"; label theta = "Long term Mean"; fit ffrate / fiml breusch=( lag_ffrate ) out=resid outresid; run;   The test for heteroscedasticity is positive and the residuals are shown in Figure 1.1. Note that the residuals are clearly heteroscedastic. To correct for the heteroscedasticity, you can add a variance model. Another popular model for this type of data is the Cox Ingersoll Ross model, which accounts for the heteroscedasticity: The following PROC MODEL statements are used to fit this model: proc model data=weeklyrates; ffrate = lag(ffrate) + kappa * (theta - lag(ffrate)); h.ffrate = sigma**2 * lag(ffrate); label kappa = "Speed of Mean Reversion"; label theta = "Long term Mean"; label sigma ="Constant part of variance"; fit ffrate / fiml out=resid outresid; run;   The residuals are shown in Figure 1.3 . The dependency on the lag of FFRATE has been removed, but the residuals remain heteroscedastic. A GARCH(1,1) model is now considered:     Figure 1.3: Residuals from Cox Ingersoll Ross Model The following PROC MODEL statements fit a GARCH(1,1) model: proc model data=weeklyrates; ffrate = lag(ffrate) + kappa * (theta - lag(ffrate)); if ( _OBS_ = 1 ) then h.ffrate = arch0 + arch1 * mse.ffrate + garch1 * mse.ffrate; else h.ffrate = arch0 + arch1 * zlag(resid.ffrate**2) + garch1 * zlag(h.ffrate); fit ffrate / fiml out=resid outresid; run;   The residuals from this estimation are shown in Figure 1.4 . They look more homoscedastic than those produced by the other models, however, some additional investigation may be necessary on the data point occurring around 1987.     Figure 1.4: Residuals from Garch(1,1)   References Chan, K.C. and Karolyi, G. Andrew and Longstaff, Francis A. and Sanders, Anthony B. (1992),"An Empirical Comparison of Alternate Models of the Short-Term Interest Rate", The Journal of Finance, 47(3), pp. 1209--1227. Greene, William H. (1993), "Econometric Analysis," New York: Macmillian Publishing Company Inc. SAS Institute Inc. (1993), SAS/ETS User's Guide, Version 8, Cary, NC: SAS Institute Inc. ... View more

Overview Heteroscedastic two-stage least squares regression is a modification of the traditional two-stage least squares used to estimate simultaneous equation models when the disturbances are heteroscedastic. One can use the MODEL procedure in SAS/ETS to compute the two-stage heteroscedastic estimates. Details Consider the general nonlinear simultaneous equations model   where  is a real valued vector function of  ,  ,  , g is the number of equations, l is the number of exogenous variables, p is the number of parameters and t ranges from 1 to n.   is a vector of instruments.  is an unobservable disturbance vector with the following properties: The nonlinear two stage least squares (N2SLS) is a commonly used single equation estimation method that consists of using instrumental variables (IV) that are uncorrelated with the disturbances to obtain predicted values for the endogenous variables. Such predicted values replace the right hand side endogenous variables in the model to obtain consistent estimates of the parameters. Formally, the N2SLS estimator θ is the value that minimizes the quantity     A key assumption of N2SLS is that the disturbances are i.i.d.. In the presence of heteroscedasticity, the 2SLS estimator is still consistent but no longer efficient. In such case, the nonlinear heteroscedastic 2SLS (NH2SLS) estimator is a single-equation estimator that attains a gain in efficency with respect to the N2SLS estimator by weighting observations with White's estimate of the error correlation matrix. Formally, the NH2SLS estimator is obtained by minimizing the following quantity: where  is a gk ×gk block-diagonal matrix with each block j given by, for j = 1, ... ,g. Note that the NH2SLS can be interpreted as a GMM estimator where the variance of the moment function V is taken to be block-diagonal. Each block, corresponding to each single equation, is estimated using the residuals from the N2SLS estimates. The following table illustrates the relationship among some commonly used estimation methods for simultaneous equation models. Table 1.1: Summary of PROC MODEL Estimation Methods Method Objective Function N2SLS NH2SLS N3SLS GMM   where m n  is first moment of the crossproduct  . Klein's Model The following is the well known Klein's (1950) model I, along with an excerpt of the data. The complete data can be found in the SAS code attached, along with the code used in this example. options ls=80; data klein; input year c p w i x wp g t k wsum; ... more datalines ... 1934 48.7 12.3 30.6 -3.0 49.7 6.0 4.0 6.8 199.0 36.6 1935 51.3 14.0 33.2 -1.3 54.4 6.1 4.4 7.2 197.7 39.3 1936 57.7 17.6 36.8 2.1 62.7 7.4 2.9 8.3 199.8 44.2 1937 58.7 17.3 41.0 2.0 65.0 6.7 4.3 6.7 201.8 47.7 1938 57.5 15.3 38.2 -1.9 60.9 7.7 5.3 7.4 199.9 45.9 1939 61.6 19.0 41.6 1.3 69.5 7.8 6.6 8.9 201.2 49.4 1940 65.0 21.1 45.0 3.3 75.7 8.0 7.4 9.6 204.5 53.0 1941 69.7 23.5 53.3 4.9 88.4 8.5 13.8 11.6 209.4 61.8 ; run; The endogenous variables are on the left hand side of the equations. The exogenous variables are government wage bills, W g , government spending, G, indirect business taxes plus net export, T, the time trend in terms of years from  1931, Yt, and the constant terms.   The following PROC MODEL statements estimate the model using 2SLS and store the estimated variance of the moment functions  in the dataset vdata via the OUTV= option.   proc model data=klein; instruments klag plag xlag wp g t yr; c=c0 + c1*p + c2*plag + c3*wsum; i=i0 + i1*p + i2*plag + i3*klag; w=w0 + w1*x + w2*xlag + w3*yr; fit c i w / 2sls outv=vdata vardef=n kernel=(bart,0,); title '2sls results'; run; quit;     The VARDEF=n and KERNEL=(bart, 0,) options have been used so that no correction for degrees of freedom or smoothing over time are performed when the matrix  is computed. See the chapter on the MODEL procedure in the SAS/ETS User's Guide for more details on these options. One can estimate the model with H2SLS using the GMM method with the variance matrix  set equal to a block diagonal matrix, where each block is given by the 2SLS estimate for the corresponding equation. See Greene (1997), p. 759 for details. The following SAS statements create the block-diagonal matrix  using the matrix stored in the dataset vdata and print the first 15 observations.   data vblkdiag; set vdata; if (eq_row ^= eq_col) then value=0; /* make V block-diagonal */ proc print data=vblkdiag(obs=15); title 'Block-diagonal V-matrix'; run;     The block-diagonal matrix  is then used in the VDATA= option of the MODEL procedure, so as to obtain the H2SLS estimates via GMM estimation. Note that all the exogenous variables are used as instruments. Also note that the NOGENGMMV option has been turned on in the FIT statement. This specifies the use of GMM variance under the optimal weighting matrix instead of the arbitrary weighting matrix. See the Estimation Methods in the SAS/ETS User's Guide. proc model data=klein; instruments klag plag xlag wp g t yr; c=c0 + c1*p + c2*plag + c3*wsum; i=i0 + i1*p + i2*plag + i3*klag; w=w0 + w1*x + w2*xlag + w3*yr; fit c i w / no2sls gmm vdata=vblkdiag vardef=n nogengmmv; title 'H2SLS results'; run; quit;   The results of the estimation are shown in Figure 1 . You can compare them with Greene (1997), p. 760. H2SLS results The MODEL Procedure Nonlinear GMM Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| c0 14.74433 1.1596 12.71 <.0001 c1 0.075792 0.0936 0.81 0.4291 c2 0.166269 0.0825 2.02 0.0599 c3 0.849365 0.0356 23.85 <.0001 i0 21.40696 6.3296 3.38 0.0035 i1 0.18586 0.1299 1.43 0.1706 i2 0.551308 0.1253 4.40 0.0004 i3 -0.16056 0.0305 -5.27 <.0001 w0 2.674615 0.7915 3.38 0.0036 w1 0.455802 0.0300 15.21 <.0001 w2 0.110765 0.0323 3.43 0.0032 w3 0.1306 0.0237 5.51 <.0001 Figure 1: H2SLS Results   References Greene, W.H. (1997), Econometric Analysis, Third Edition, Upple Saddle River, NJ: Prentice Hall, Inc. Klein, L. (1950), Economic Fluctuations in the Unites States 1921-1940. New York: John Wileys and Sons. SAS Institute, Inc. (1999), SAS/ETS User's guide, Version 8, Volume 2, Cary, NC: SAS Institute, Inc. SAS Institute, Inc. (2003), SAS/ETS User's guide, Version 9, Cary, NC: SAS Institute, Inc. ... View more

Overview Figure 1: Comparison of Density Estimates for the ARCH 1 Parameter Estimator Under Different Imputation Methods Missing data in time series constitutes a problem that is not easily solved. Standard techniques available for cross-sectional data may not be applicable due to the dynamic nature of the models. Ad hoc procedures may either ignore some important facts or be difficult to implement. Multiple imputation has been used for several years with cross sectional models to deal with missing data. Its flexibility makes it an adequate tool for times series data as well. In this example, you will see how to implement a simple form of multiple imputation for time series to fit a GARCH(1,1) model when some of the data are missing. The only tools that you will need are the MODEL procedure, the MIANALYZE procedure, and some DATA Step statements.   Details Missing observations in times series data can be a difficult problem to overcome. Due to the dynamic nature of the models, simple techniques in use for cross-sectional models, like row wise deletion, are not applicable or may lead to questionable results. One common way to deal with it is to employ a Kalman filter (Kohn and Ansley 1986). However, this usually requires the additional task of rewriting the model in state-space form, which may be difficult or may not be supported by the currently available software, if, for example, the model is nonlinear, like a GARCH model. Single imputation (SI) is another common way to deal with missing data. Each missing value is replaced by an imputed value, for example by interpolation. The replaced values are then treated as if they were actual data values. This approach enables analysis with procedures designed for complete data sets. While this method is simple and can be applied to virtually every data set, it does not account for the uncertainty about the predictions of the imputed values. Therefore, the estimated variances of the parameters are biased toward zero, leading to statistically invalid inferences, in the sense of Rubin (1987). Multiple Imputation (MI) is a methodology for dealing with incomplete data that has received considerable attention in the last couple of decades (see, for example, Rubin 1987 and 1996, and Schafer 1997.) It maintains the flexibility and relative ease of application of single imputation, while taking into account the variability due to the imputation of the missing values. Instead of filling each missing value with only one single value, MI imputes a set of m plausible values that represent the uncertainty about the right value. This approach creates m imputed replications of the original data sets that can be analyzed by standard procedures for complete data sets. If Q is the quantity of interest, then the analysis of the m imputed data sets will produce m estimates  of Q, with the corresponding estimates  of their variance. The overall estimate  of Q is the average of the m complete-data estimates  : As shown by Rubin (1987), the variance V of the multiple imputation point estimate is the average of the estimated variances within each complete data set,  , plus the sample variance of the point estimates  between the datasets, multiplied by a correction factor. where is the sample variance of the estimated  between the datasets. You can refer to "The MI Procedure" and "The MIANALYZE Procedure" chapters in the SAS/STAT User's Guide for more details on multiple imputation. The multiple imputation literature so far has focused mainly on applications in cross-sectional models for survey data. Only recently, panel and time series data have been considered (see King et al. 2001, and Hopke, Chuanhai, and Rubin 2001.) In the same fashion, software for multiple imputation, like the MI procedure, supports mainly cross-sectional models. Although, from a theoretical point of view, there is no reason not to use multiple imputation for time series data, its application has been difficult in practice. In the following example you will see how to use the Monte Carlo feature within the MODEL procedure to implement a simple form of multiple imputation, and to estimate a zero mean GARCH(1,1) model with missing data. Other algorithms of multiple imputation (see Hopke, Chuanhai, and Rubin 2001) may be more appropriate, however, this method is of very simple implementation, and provides results that are statistically superior to simple imputation in the case that is presented. The MIANALYZE procedure can then combine the results from the different imputed data. The procedure has been designed to provide for great flexibility and enables you to analyze output obtained using a BY statement from virtually any SAS procedure. It reads parameter estimates and associated standard errors, and derives valid inferences according to the preceding formulas. The GARCH(1,1) model with mean zero is described by the following equations:   The parameter vector θ is given by  . For this example, a total of 1,000 data points are simulated from the above model, where the true values of the parameters are  ,  , and  . Then, about 10% of the original data are randomly set to missing. (The full code can be found in the attached SAS source file.) The estimation procedure can be divided into four main steps. The first step is to find some plausible preliminary estimates of the model parameters. One way to accomplish this is to use single imputation by replacing the missing observation with the mean, which, in this case, is zero. Other possible options include single imputation with polynomial interpolation, or using the longest set of data that does not contain missing observations.   data single; set missing; if y=. then y=0; /* replace y with mean */ run; The following PROC MODEL statements estimate the parameters of the GARCH(1,1) model using the dataset in which the missing data are replaced by the expected value. The estimated parameter vector  is displayed in Figure 2 .     proc model data=single; parms arch0 arch1 garch1; /* Mean zero model --------*/ y = 0; /* Garch Variance ---------*/ h.y = arch0 + arch1 * xlag( resid.y**2, mse.y) + garch1 * xlag( h.y, mse.y ) ; /* Fit the garch model ----*/ fit y / data=single method=marquardt fiml outest=prelest outcov outs=s prl=both; bounds arch0 arch1 garch1 >= 0; title 'Single imputation'; run; Single imputation The MODEL Procedure Nonlinear FIML Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| arch0 0.60608 0.1619 3.74 0.0002 arch1 0.24358 0.0458 5.31 <.0001 garch1 0.472489 0.1003 4.71 <.0001 Parameter Wald 95% Confidence Intervals Parameter Value Lower Upper arch0 0.6061 0.2888 0.9233 arch1 0.2436 0.1537 0.3334 garch1 0.4725 0.2759 0.6691 Parameter Likelihood Ratio 95% Confidence Intervals Parameter Value Lower Upper arch0 0.6061 0.4855 0.9233 arch1 0.2436 0.2094 0.3334 garch1 0.4725 0.3978 0.5472 Figure 2: Single Imputation Estimates The estimate  , and its estimated covariance matrix, are saved in the prelest data set. Likewise, the estimated standard error of the residual is saved in the s data set. The second step is to generate the imputed data. First define how many you want. In this case, 20 sets of imputations are requested.   %let m=20; /* number of imputations */ Then generate them with the Monte Carlo feature within the MODEL procedure. The following SOLVE statement creates m replications of the input data with missing values replaced by imputed data. It accomplishes this by drawing m vectors  of parameters from a multivariate normal distribution whose mean and covariance matrix are the estimated parameter vector  and covariance matrix  . Bollerslev (1986) proved that the asymptotic distribution of the maximum likelihood estimators of a GARCH model is multivariate normal. Therefore, this choice of distribution seems appropriate. The errors et are simulated from a univariate normal distribution with mean zero, while the standard error is taken from the s data set. The SOLVE statement supports simulation from several other distributions. Refer to the "Details: Simulation" section of "The Model Procedure" chapter in the SAS/ETS User's Guide for more details.   /* Impute the missing point in DATA=missing data set */ solve y / data=missing estdata=prelest sdata=s random=&m seed=123 out=monte forecast;   The 20 sets of imputed data values are saved in the monte output data set and are identified by the _REP_ variable. The data corresponding to _REP_=0 are the original data values. The following output in Figure 3  shows that the last two missing values on the original data set are replaced by imputed values in the first two imputed replications of the data set, while the rest of the nonmissing data is kept unchanged. That is accomplished by the FORECAST option of the SOLVE statement.   Original data set Obs y 989 2.18025 990 -1.97288 991 . 992 -1.29694 993 -1.48274   First set of imputed data Obs _REP_ y 1989 1 2.18025 1990 1 -1.97288 1991 1 -3.74653 1992 1 -1.29694 1993 1 -1.48274   Second set of imputed data Obs _REP_ y 2989 2 2.18025 2990 2 -1.97288 2991 2 -2.94843 2992 2 -1.29694 2993 2 -1.48274 Figure 3: Multiply Imputed Data Sets The next step is to reestimate the model for each set of imputed values using the BY statement. The starting values of the parameters are the estimates of the single imputation. The original set of data, corresponding to _REP_=0, is dropped from the input data set. It is important to note that the estimated model is the same as the one used for imputation. In the words of Shafer (1997), the imputer's model coincides with the analyst's model, so that any problem arising from a discrepancy is avoided. Refer to "The MI Procedure" chapter in the SAS/STAT User's Guide and Schafer (1997) for an in-depth discussion about the choice of the imputer's model versus the analyst's model.   /* Fit again y with the imputed data */ fit y / data=monte(where=(_rep_>0)) method=marquardt fiml estdata=prelest sdata=s nools outest=imp_est outcov; bounds arch0 arch1 garch1 >= 0; by _rep_; title 'Multiple Imputation'; run;   The NOOLS option skips the preliminary OLS estimation that PROC MODEL does by default to find starting values for the parameters, since here the starting values can be provided by the single imputation estimates. The last step is to summarize the results with the MIANALYZE procedure. The OUTEST= data set of PROC MODEL needs to be converted into a type=EST data set, and the name of the BY variable needs to be _Imputation_ in order to be processed by MIANALYZE.   data mi_in_est(type=EST); set imp_est; rename _rep_ = _Imputation_; if _name_="" then _type_="PARMS"; else _type_="COV"; run; The following MIANALYZE statements produce parameter estimates and test statistics for the parameters of interest. The relevant part of the output is shown in Figure 4. The EDF= option specifies the complete-data degrees of freedom. Refer to "The MIANALYZE Procedure" chapter in the SAS/STAT User's Guide for more details on this option.     proc mianalyze data=mi_in_est edf=1000; modeleffects arch0 arch1 garch1; run;     The MIANALYZE Procedure Multiple Imputation Parameter Estimates Parameter Estimate Std Error Std Error 95% Confidence Limits DF arch0 0.680916 0.198299 0.198299 0.289496 1.072336 171.59 arch1 0.301525 0.061156 0.061156 0.181066 0.421983 245.08 garch1 0.422048 0.117238 0.117238 0.190272 0.653823 140.69 Figure 4: Multiple Imputation Parameter Estimates To verify that this method produces results that are statistically superior to single imputation, a Monte Carlo experiment was conducted. In this experiment, 10,000 time series were simulated from a GARCH(1,1) model, and the nominal 95% confidence intervals were computed both in the case of single imputation where the missing values were replaced either by the expected value or by polynomial interpolation, and in the case of multiple imputation according to the scheme outlined above. The results are shown in Figure 5 . The type variable indicates the type of confidence intervals: Wald for the complete data ('wald_full'), likelihood ratio for the complete data ('likl_full'), Wald computed with single imputation by mean replacement ('wald_mean'), likelihood ratio computed with single imputation by mean replacement ('likl_mean'), Wald computed with single imputation by interpolation ('wald_inte'), likelihood ratio computed with single imputation by interpolation ('likl_inte'), or with multiple imputation ('mi'). The actual coverage of the confidence intervals obtained with multiple imputation is much closer to the nominal coverage than the one of confidence intervals obtained with single imputation. Therefore, multiple imputation according the the procedure just outlined leads to better statistical inference, in the sense of Rubin (1987), for this particular model. Figure 1  shows also that the distribution of mi estimates of the arch1 parameter is closer in shape and location to the distribution of the estimates when the complete data set is used. Parameter=arch0   Obs type coverage 1 wald_full 0.9406 2 likl_full 0.7236 3 mi 0.9445 4 wald_mean 0.8501 5 likl_mean 0.7245 6 wald_inte 0.7458 7 likl_inte 0.5687 Parameter=arch1   Obs type coverage 8 wald_full 0.9452 9 likl_full 0.7512 10 mi 0.9136 11 wald_mean 0.7652 12 likl_mean 0.7186 13 wald_inte 0.7412 14 likl_inte 0.6423 Parameter=garch1   Obs type coverage 15 wald_full 0.9359 16 likl_full 0.8173 17 mi 0.9409 18 wald_mean 0.8616 19 likl_mean 0.7789 20 wald_inte 0.7554 21 likl_inte 0.6635 Figure 5: Actual Coverage of 95% Confidence Intervals The source code for the Monte Carlo study is available upon request. The %Distribute macro system of Doninger and Tobias (2001) has provided the tools to perform the study within a reasonable amount of time.   References Bollerslev, T. (1986), "Generalized Autoregressive Conditional Heteroskedasticity," Journal of Econometrics, 31, 307--327. Doninger, C. and Tobias, R. (2001), "The %Distribute System for Large-Scale Parallel Computation in the SAS System.". Cary, NC: SAS Institute, Inc. Accessed 12 October 2001. Hopke, P. K., Chuanhai, L., and Rubin, D. B. (2001), "Multiple Imputation for Multivariate Data with Missing And Below-Treshold Measurement: Time-Series Concentration of Pollutants in the Artic," Biometrics, 57, 22--33. King, G., Honaker, J., Joseph, A., and Scheve, K. (2001), "Analyzing Incomplete Political Science Data: An Alternative Algorithm for Multiple Imputation," American Political Science Review, 95, 49--69. Kohn, R. and Ansley, C.F (1986), "Efficient Estimation and Prediction in Time Series Regression Models," Biometrica, 72, 694--697. Rubin, D.B. (1987), Multiple Imputation for Nonresponse in Surveys, New York: John Wiley & Sons, Inc. Rubin, D. (1996), "Multiple Imputation after 18+ Years," Journal of the American Statistical Association, 91, 473--489. Schafer, J.L. (1997), Analysis of Incomplete Data, London: Chapman and Hall. SAS Institute, Inc. (2000), SAS/ETS User's Guide, Version 8.2, Cary, NC: SAS Institute, Inc. SAS Institute, Inc. (2000), SAS/STAT User's Guide, Version 8.2, Cary, NC: SAS Institute, Inc. ... View more

Overview It is said that forecasting is much more of an art than a science. Many practitioners use a visual comparison of several forecasts to assess their accuracy and choose among the various forecasting methods. This example shows how to plot on the same graph the original values of a time series variable and the predicted values from three different forecasting methods, thus facilitating a visual comparison. Analysis This example uses lead production data as the forecast variable. The values represent monthly totals of U.S. lead production, in tons, for the period January 1990 to September 1992. Begin by resetting all graphics options and entering the lead production data with a SAS DATA step.   goptions reset=all; /* ------- Read Initial Data ------- */ data leadprd; input date:monyy5. leadprod @@; format date monyy5.; title 'Lead Production Data'; title2 '(in tons)'; datalines; jan90 38500 feb90 37900 mar90 36900 apr90 38600 may90 36400 jun90 33300 jul90 34000 aug90 38000 sep90 37400 oct90 42300 nov90 36900 dec90 34800 jan91 33900 feb91 34000 mar91 37200 apr91 33300 may91 29800 jun91 24700 jul91 30800 aug91 31100 sep91 32400 oct91 32900 nov91 29100 dec91 31800 jan92 32100 feb92 30500 mar92 36800 apr92 30300 may92 29500 jun92 24700 jul92 27600 aug92 23800 sep92 21400 ; run;     Next produce your forecasts and save their predicted values to SAS data sets. This example uses the forecasting capabilities of the FORECAST, the ARIMA, and the REG procedures. The OUT1STEP option of PROC FORECAST specifies that only the one-step-ahead forecasts are output to the data set LEADOUT1. The LEAD= option produces forecasts for 12 months beyond the sample period.     /* -- Using PROC FORECAST -- */ proc forecast data=leadprd out=leadout1 out1step lead=12 interval=month; id date; var leadprod; run;     The F statement in PROC ARIMA is used to generate forecasts based on the ARIMA model estimated with the E statement. The LEAD= option, as in the FORECAST procedure, produces forecasts for 12 months beyond the sample period.   /* -- Using PROC ARIMA -- */ proc arima data=leadprd; i var=leadprod nlag=15; e p=1; f lead=12 interval=month id=date out=leadout2; run; quit;     To estimate a time trend for the lead prediction data, it is necessary to create a new variable T that spans both the sample and forecast periods.   /* -- Using PROC REG (Time Trend Regression) -- */ data ttrend; set leadout2; t+1; run;     Use PROC REG to create an output data set, LEADOUT3, containing the predicted values from a regression of LEADPROD on a time trend.   proc reg data=ttrend; model leadprod = t; output out=leadout3 p=ptrend; run; quit;     The key to overlaying plots of time series forecasts is in the data management. Use a DATA step to collect all of the forecasts in a data set named FINAL. Use the KEEP= and the RENAME= options to combine the three output data sets.   /* ------- Data Management ------- */ data final; set leadout1(keep=date leadprod rename=(leadprod=pfore)); set leadout2(keep=date leadprod forecast rename=(leadprod=actual forecast=parima)); set leadout3(keep=date ptrend); run;     Use graphics options along with SYMBOL, AXIS, and LEGEND statements to provide a more polished look to your graph.   /* ------- Graphics Output ------- */ goptions cback=white ctitle=bl ctext=bl border ftitle=centx ftext=centx; title 'Lead Production Data'; title2 'Plot of Forecast for Lead Production'; symbol1 i=spline width=1 v=dot c=black; /* for actual */ symbol2 i=spline width=2 v=none c=red; /* for pfore forecast */ symbol3 i=spline width=2 v=none c=green; /* for parima forecast */ symbol4 i=spline width=2 v=none c=blue; /* for ptrend forecast */ axis1 offset=(1 cm) label=('Year') minor=none order=('01jan90'd to '01jan94'd by year); axis2 label=(angle=90 'Lead Production') order=(20000 to 45000 by 5000); legend1 across=1 cborder=black position=(top inside right) offset=(-2,0) value=(tick=1 'ACTUAL' tick=2 'PROC FORECAST' tick=3 'PROC ARIMA' tick=4 'TIME TREND') shape=symbol(2,.25) mode=share label=none;     Use the GPLOT procedure with the OVERLAY option to plot the actual time series data along with the predictions from the three forecasting methods. The HREF= option creates a reference line perpendicular to the horizonal axis at September 1, 1992, the dividing date between the sample and forecast periods.   proc gplot data=final; format date year4.; plot actual * date = 1 pfore * date = 2 parima * date = 3 ptrend * date = 4 / overlay noframe href='01sep92'd chref=red vaxis=axis2 vminor=1 haxis=axis1 legend=legend1; run; quit;     The graph produced is shown in Figure 1    Figure 1: Forecasts of Lead Production   References SAS Institute Inc. (1996), Forecasting Examples for Business and Economics Using the SAS System, Cary, NC: SAS Institute Inc. SAS Institute Inc. (1993), SAS/ETS User's Guide, Version 6, Second Edition, Cary, NC: SAS Institute Inc. SAS Institute Inc. (1990), SAS/GRAPH Software: Reference, Version 6, First Edition, Volume 2, Cary, NC: SAS Institute Inc. ... View more

Overview Graphic visualization of time series variables is helpful in identifying and interpreting relationships in data. High resolution color graphics plots convey the information quickly and concisely. This example makes use of the ODS graphics for PROC ARIMA and compares them to the graphics produced using the GPLOT procedure. The user will be able to plot forecasts and the associated 95% confidence band for an intervention model of driver casualties (deaths and serious injuries) in road accidents in Great Britain. The intervention variable is the passage of a mandatory seatbelt law in January, 1983. Analysis In January 1983, the British government passed a mandatory seatbelt law in order to reduce the number of casualties in road accidents. Seatbelt laws are often controversial. Opponents sometimes claim that they do not significantly reduce the number of casualties. A graphical analysis of the forecasted values and the upper and lower 95% confidence limits from an intervention model gives us valuable insights into the impact of the legislation. The driver casualties are taken from Harvey (1989). The INJURIES variable contains the number of deaths and serious injuries to car drivers on roads in Great Britain from January 1980 to December 1984. Twelve missing values are added to the end of the data set for use in the example. In January 1983, a mandatory seatbelt law took effect. The BELTLAW variable takes on the value 0 before January 1, 1983, and a value of 1 thereafter.   data drivers; format date monyy.; input date:monyy5. injuries @@; beltlaw = (date ge '01jan83'd); datalines; jan80 1665 feb80 1361 mar80 1506 apr80 1360 may80 1453 jun80 1522 jul80 1460 aug80 1552 sep80 1548 oct80 1827 nov80 1737 dec80 1941 jan81 1474 feb81 1458 mar81 1542 apr81 1404 ... more lines ... sep84 1444 oct84 1575 nov84 1737 dec84 1763 jan85 . feb85 . mar85 . apr85 . may85 . jun85 . jul85 . aug85 . sep85 . oct85 . nov85 . dec85 . ; run;   The intervention model is estimated with the ARIMA procedure. Briefly, the procedure uses the input data set DRIVERS. The I statement identifies the variable INJURIES, the E statement estimates the ARIMA model specified, and the F statement outputs monthly forecasts to the output data set, FORE2. For a more thorough analysis of the intricacies of working with ARIMA models, see forecasting a Seasonal ARMA Process. Before the introduction of ODS graphics, the GPLOT procedure had been used on many different levels. Simple plots can be produced with ease; for instance, the time series plot in Figure 1.1  of the variable INJURIES can be produced with the following statement:   proc gplot data=drivers; plot injuries*date=1; symbol1 v=star c=blue; title "Time Series Plot"; run; quit; title;   Figure 1.1 Simple Time Series Plot of INJURIES   In the past, to achieve professional quality graphics, you would have run the ARIMA procedure, saved the output into a DATA set, and used the output in the GPLOT procedure. This would have required some formatting to achieve professional quality graphics. The following code produces the forecast plot in Figure 1.2  using the ARIMA and GPLOT procedures. proc arima data=drivers; i var=injuries crosscorr=beltlaw noprint; e input=beltlaw p=(1)(2 12); f lead=12 out=fore2 id=date interval=month; run; quit; goptions cback=white ctitle=bl ctext=bl ftitle=centx ftext=centx reset=(axis symbol); proc gplot data=fore2; format date year4.; plot injuries*date=1 forecast*date=2 l95*date =3 u95*date =3 / overlay haxis=axis1 vaxis=axis2 vminor=4 href='01jan83'd lh=2; title 'Driver Casualties Data'; title2 'Monthly Totals'; axis1 offset=(1 cm) label=('Year') order=('01jan80'd to '01jan86'd by year); axis2 label=(angle=90 'Casualties') order=(750 to 2250 by 500); symbol1 i=join l=1 c=red; symbol2 h=2 pct v=star c=blue; symbol3 l=20 i=join c=green; footnote1 c=r f=centx ' --- Actual' c=b f=centx ' * Forecast' c=g f=centx ' --- U95/L95'; run; quit; title; footnote; Figure 1.2 Forecast Plot Produced by PROC GPLOT   It is always good form to specify the formating options for your graph in the GOPTIONS statement. The CBACK= option sets the color of the background, while the CTITLE= and CTEXT= options enable you to choose a font from among the many supported by SAS/GRAPH software for the title and text. In this example, the CENTX (Century Expanded) font is used for both the title and text. The RESET= option resets the title and the text. The axes can be formatted using the AXISn definitions. The OFFSET= option specifies the amount of space to offset the first and last major tick marks from the origin and end of the axis. The LABEL= option enables you to label each axis in the graph. The ANGLE=90 option in the AXIS2 statement rotates the label 90 degrees. The ORDER= option specifies the data values in the order in which they are to appear on the axis. The SYMBOL definitions are used to specify the shape (V= option), as well as the size (H= ) and color (C= ), of the plot symbols used to mark the data points. The I= option enables you to interpolate between data points by a number of methods. The type of line used is specified with the L= option. A value of L=1 corresponds to a solid line, while L=20 is a dashed line. Footnotes can be used to provide additional information. In this example, a footnote is used to add a legend to the plot. The GPLOT procedure uses the input data set FORE2. The OVERLAY option in the PLOT statement plots the time series INJURIES, FORECAST, L95, and U95 on the same graph using the symbols indicated. The HAXIS and VAXIS options specify the horizontal and vertical axes to be used in the plot. The VMINOR=4 option places four minor tick marks between each major tick mark on the vertical axis. The HREF= option draws a horizontal reference line at January 1, 1986, the date the mandatory seatbelt law took effect, and the LH= option specifies the line type. With the introduction of ODS Graphics, the user need only turn on the graphics feature to have several standard graphics displayed. The statement ODS GRAPHICS ON; needs to be specified to enable this feature. In Table 1.1 , you can see an overview of the different graphs that are produced from the different ARIMA statements. More information about the ODS graphics in PROC ARIMA can be found in the SAS/ETS User's Guide, Version 9.2.   Table 1.1 ODS Graphics Produced by PROC ARIMA Associated Statement Plot Description FORECAST Forecast Plot ESTIMATE Residual Autocorrelation Plot ESTIMATE Residual Inverse Autocorrelation Plot ESTIMATE Residual Partial Autocorrelation Plot ESTIMATE Residual White Noise Probability Plot ESTIMATE Residual Histogram ESTIMATE Residual Normal Q-Q Plot ESTIMATE Residual Partial Autocorrelation Plot IDENTIFY Series Autocorrelation Plot IDENTIFY Series Inverse Autocorrelation Plot IDENTIFY Series Partial Autocorrelation Plot IDENTIFY Trend Analysis Plot IDENTIFY Cross Correlation Plot   For better control of graphics, you can also specify global plot option ONLY in the PLOTS= option in the PROC ARIMA statement. This suppresses generation of default plots and produces plots specified in the subsequent FORECAST or any other plot options. For example, specific plot requests are given by specifying a statement as follows: PLOTS(ONLY)=(SERIES(CORR CROSSCORR) RESIDUAL(NORMAL SMOOTH)) The SERIES(CORR CROSSCORR) option produces a panel of plots that are useful for the trend and correlation analysis of the series and cross-correlation plots respectively. The RESIDUAL(NORMAL SMOOTH) option produces the histogram and Q-Q plot with the NORMAL option and a scatter plot overlaid with a smooth fit with the SMOOTH option. You can suppress the generation of all plots by using the PLOTS=NONE option. To produce all plots relevant for a particular analysis, use the PLOTS=ALL option. For more details about the various options available with ODS graphics, see SAS/ETS User's Guide, Version 9.2. You will need to turn on the ODS graphics as shown below. The following code produces the ODS Forecast Plot in Figure 1.3. The PLOTS(ONLY)= FORECAST(FORECAST) option produces a plot that shows the one-step ahead forecasts as well as the multi-step ahead forecasts. ods graphics on; proc arima data=drivers plots(only)=(forecast(forecast)); i var=injuries crosscorr=beltlaw; e input=beltlaw p=(1)(2 12); f lead=12 out=fore2 id=date interval=month; run; quit; Figure 1.3 Forecast Plot Produced with ODS Graphics   From the graphs, it appears that the passage of the seatbelt law caused a significant decrease in the number of road casualties. References Harvey, A.C. (1989), Forecasting, Structural Time Series Models and the Kalman Filter, Cambridge: Cambridge University Press. SAS Institute Inc. (1996), Forecasting Examples for Business and Economics Using the SAS System, Cary, NC: SAS Institute Inc. SAS Institute Inc. (1993), SAS/ETS User’s Guide, Version 6, Second Edition, Cary, NC: SAS Institute Inc. SAS Institute Inc. (2008), SAS/ETS User’s Guide, Version 9.2, Second Edition, Cary, NC: SAS Institute Inc. SAS Institute Inc. (1990), SAS/GRAPH Software: Reference, Version 6, First Edition, Volume 2, Cary, NC: SAS Institute Inc. SAS Institute Inc. (2004), SAS/GRAPH Software: Reference, Version 9, Cary, NC: SAS Institute Inc. ... View more

Overview One of the classical assumptions of the ordinary regression model is that the disturbance variance is constant, or homogeneous, across observations. If this assumption is violated, the errors are said to be "heteroscedastic." Heteroscedasticity often arises in the analysis of cross-sectional data. For example, in analyzing public school spending, certain states may have greater variation in expenditure than others. If heteroscedasticity is present and a regression of spending on per capita income by state and its square is computed, the parameter estimates are still consistent but they are no longer efficient. Thus, inferences from the standard errors are likely to be misleading. Testing for Heteroscedasticity There are several methods of testing for the presence of heteroscedasticity. The most commonly used is the Time-Honored Method of Inspection (THMI). This test involves looking for patterns in a plot of the residuals from a regression. Two more formal tests are White's General test (White 1980) and the Breusch-Pagan test (Breusch and Pagan 1979). The White test is computed by finding nR 2 from a regression of e i 2 on all of the distinct variables in  , where X is the vector of dependent variables including a constant. This statistic is asymptotically distributed as chi-square with k-1 degrees of freedom, where k is the number of regressors, excluding the constant term. The Breusch-Pagan test is a Lagrange multiplier test of the hypothesis that the independent variables have no explanatory power on the e i 2 's. If u equals (e 1 2 ,e 2 2 , . . . ,e n 2 ), i equals an n ×1 column of ones, and  , then Koenkar and Bassett's (1982) robust variance estimator   computes the test statistic as which is asymptotically distributed as chi-square with degrees of freedom equal to the number of variables in Z.   Correcting for Heteroscedasticity One way to correct for heteroscedasticity is to compute the weighted least squares (WLS) estimator using an hypothesized specification for the variance. Often this specification is one of the regressors or its square. This example uses the MODEL procedure to perform the preceding tests and the WLS correction in an investigation of public school spending in the United States.   Analysis If y is public school spending and x is per capita income, and assuming that the variance of the error term is proportional to xi2, then the regression model in this example can be written as where i = 1, ... ,51 is a state index. The sample consists of 51 observations of per capita expenditure on public schools and per capita income for each state and the District of Columbia in 1979. The following DATA step reads in the 51 observations, transforms the variable INC by multiplying it by 10-4 (for consistency with Greene 1993), creates the variable INC2 as the square of income, and then deletes Wisconsin from the sample due to a missing value for expenditure.   data hetero1; input st exp inc; inc=inc/10000; inc2=inc**2; if exp = . then delete; datalines; 1 275 6247 2 275 6183 3 531 8914 ... ; run;   Testing for Heteroscedasticity You can use the MODEL procedure for the initial investigation of the model. The following commands estimate the preceding model, perform two different tests for heteroscedasticity (the White and the Breusch-Pagan), and output the residuals into a data set for further investigation.   proc model data=hetero1; parms a1 b1 b2; exp = a1 + b1 * inc + b2 * inc2; fit exp / white pagan=(1 inc inc2) out=resid1 outresid; run; quit;     Ordinary Least Squares   The MODEL Procedure Nonlinear OLS Summary of Residual Errors  Equation DF Model DF Error SSE MSE Root MSE R-Square Adj R-Sq exp 3 47 150986 3212.5 56.6785 0.6553 0.6407   Nonlinear OLS Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| a1 832.9144 327.3 2.54 0.0143 b1 -1834.2 829.0 -2.21 0.0318 b2 1587.042 519.1 3.06 0.0037   Number of Observations Statistics for System Used 50 Objective 3020 Missing 0 Objective*N 150986 Heteroscedasticity Test Equation Test Statistic DF Pr > ChiSq Variables exp White's Test 21.16 4 0.0003 Cross of all vars   Breusch-Pagan 15.83 2 0.0004 1, inc, inc2   The estimates for the constant term and the coefficients of INC and INC2 and their associated p-values are 832.91 (0.014), -1834.20 (0.032), and 1587.04 (0.004), respectively, which all appear to be different from 0 at generally accepted levels of statistical significance. Notice, however, that both the White test (21.16) and the Breusch-Pagan test (15.83) reject the null hypothesis of no heteroscedasticity. This implies that the standard errors of the parameter estimates are incorrect and, thus, any inferences derived from them may be misleading. A plot of the residuals shows more variance in the errors of higher income states.   Correcting for Heteroscedasticity If the form of the variance is known, the WEIGHT= option can be specified in the MODEL procedure to correct for heteroscedasticity using weighted least squares (WLS). The following statement performs WLS using 1/(INC2) as the weight.   proc model data=hetero1; parms a1 b1 b2; inc2_inv = 1/inc2; exp = a1 + b1 * inc + b2 * inc2; fit exp / white pagan=(1 inc inc2); weight inc2_inv; run; quit;   Weighted Least Squares   The MODEL Procedure Nonlinear OLS Summary of Residual Errors  Equation DF Model DF Error SSE MSE Root MSE R-Square Adj R-Sq exp 3 47 238308 5070.4 71.2067 0.5983 0.5812   Nonlinear OLS Parameter Estimates Parameter Estimate Approx Std Err t Value Approx Pr > |t| a1 664.5845 333.6 1.99 0.0522 b1 -1399.28 872.1 -1.60 0.1153 b2 1311.345 563.7 2.33 0.0244   Number of Observations Statistics for System Used 50 Objective 4766 Missing 0 Objective*N 238308 Sum of Weights 91.0533       Heteroscedasticity Test Equation Test Statistic DF Pr > ChiSq Variables exp White's Test 9.31 4 0.0538 Cross of all vars   Breusch-Pagan 5.23 2 0.0733 1, inc, inc2     The corrected estimates for the constant term and the coefficients of INC and INC2 and their associated p-values are 664.58 (0.052), -1399.28 (0.115), and 1311.35 (0.024), respectively. The significance of the estimates is greatly reduced, obscuring the individual effects of the explanatory variables. The White test (9.31) and the Breusch-Pagan test (5.23) are no longer significant at the 5% level. All of the preceding calculations can be found in Greene (1993, chapter 14).   References Breusch, T. and Pagan, A. (1979), ``A Simple Test for Heteroscedasticity and Random Coefficient Variation," Econometrica, 47, 1287-1294. Greene, W.H. (1993), Econometric Analysis, Second Edition, New York: Macmillan Publishing Company. Koenkar, R., and Basset, G. (1982), ``Robust Tests for Heteroscedasticity Based on Regression Quantiles," Econometrica, 50, 43-61. SAS Institute Inc. (1993), SAS/ETS User's Guide, Version 6, Second Edition, Cary, NC: SAS Institute Inc. White, H. (1980), ``A Heteroscedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroscedasticity," Econometrica, 48, 817-838.   ... View more