Overview   This example uses the MCMC procedure to fit a Bayesian linear regression model with standardized covariates. It shows how the random walk Metropolis sampling algorithm struggles when the scales of the regression parameters are vastly different. It also illustrates that the sampling algorithm performs quite well when the covariates are standardized. Analysis   The following data set contains salary and performance information for Major League Baseball players (excluding pitchers) who played at least one game in both the 1986 and 1987 seasons. The salaries are for the 1987 season (Time Inc.; 1987) , and the performance measures are from the 1986 season (Reichler; 1987) . data baseball; input logSalary no_hits no_runs no_rbi no_bb yr_major cr_hits @@; yr_major2 = yr_major*yr_major; cr_hits2 = cr_hits*cr_hits; label no_hits="Hits in 1986" no_runs="Runs in 1986" no_rbi="RBIs in 1986" no_bb="Walks in 1986" yr_major="Years in MLB" cr_hits="Career Hits" yr_major2="Years in MLB^2" cr_hits2="Career Hits^2" logSalary = "log10(Salary)"; datalines; . 66 30 29 14 1 66 2.6766936096 81 24 38 39 14 835 ... more lines ... 2.84509804 127 65 48 37 5 806 2.942008053 136 76 50 94 12 1511 2.5854607295 126 61 43 52 6 433 2.982271233 144 85 60 78 8 857 3 170 77 44 31 11 1457 ;   The MEANS procedure produces summary statistics for these data. Summary measures are saved to the SUM_BASEBALL data set for future analysis. proc means data = baseball mean stddev; output out=sum_baseball(drop=_type_ _freq_); run;   Figure 1 displays the results.   Figure 1 PROC MEANS Summary   The MEANS Procedure   Variable Label Mean Std Dev logSalary no_hits no_runs no_rbi no_bb yr_major cr_hits yr_major2 cr_hits2 log10(Salary) Hits in 1986 Runs in 1986 RBIs in 1986 Walks in 1986 Years in MLB Career Hits Years in MLB^2 Career Hits^2 2.5730111 103.5576324 52.3333333 49.4485981 39.8878505 7.6292835 736.7570093 82.0716511 933094.76 0.3864602 44.1548123 25.0098314 25.5044973 21.1216686 4.8928790 625.7024129 93.8924809 1371125.94     Bayesian Linear Regression Model   Suppose you want to fit a Bayesian linear regression model for the logarithm of a player’s salary with density as follows:   ( 1 )     where   is the vector of covariates listed as   for i = 1,..., n = 361 baseball players. Pete Rose was an extreme outlier in 1986, and his information greatly skews results. He is omitted from this data set and analysis. The likelihood function for the logarithm of salary and the corresponding covariates is ( 2 ) where   denotes a conditional probability density. The normal density is evaluated at the specified value of log(SALARY i ) and the corresponding mean parameter μ i  defined in Equation 1. The regression parameters in the likelihood are β 0 through β 8 . Suppose the following prior distributions are placed on the parameters:     ( 3 )     where   indicates a prior distribution and   is the density function for the inverse-gamma distribution. Priors of this type with large variances are often called diffuse priors. Using Bayes’ theorem, the likelihood function and prior distributions determine the posterior distribution of the parameters as follows:     PROC MCMC obtains samples from the desired posterior distribution. You do not need to specify the exact form of the posterior distribution. The following SAS statements use the likelihood function and prior distributions to fit the Bayesian linear regression model. The PROC MCMC statement invokes the procedure and specifies the input data set. The NBI= option specifies the number of burn-in iterations. The NMC= option specifies the number of posterior simulation iterations. The SEED= option specifies a seed for the random number generator (the seed guarantees the reproducibility of the random stream). The PROPCOV=QUANEW option uses the estimated inverse Hessian matrix as the initial proposal covariance matrix.   ods graphics on; proc mcmc data=baseball nbi=50000 nmc=10000 seed=1181 propcov=quanew; array beta[9] beta0-beta8; array data[9] 1 no_hits no_runs no_rbi no_bb yr_major cr_hits yr_major2 cr_hits2; parms beta: 0; parms sig2 1; prior beta: ~ normal(0,var = 1000); prior sig2 ~ igamma(shape = 3/10, scale = 10/3); call mult(beta, data, mu); model logsalary ~ n(mu, var = sig2); run; ods graphics off;     Each of the two ARRAY statements associates a name with a list of variables and constants. The first ARRAY statement specifies names for the regression coefficients. The second ARRAY statement contains all of the covariates. The first PARMS statement places all regression parameters in a single block and assigns them an initial value of 0. The second PARMS statement places the variance parameter in a separate block and assigns it an initial value of 1. The first PRIOR statement assigns the normal prior to each of the regression parameters. The second PRIOR statement assigns the inverse-gamma prior distribution to σ 2 . The CALL statement uses the MULT matrix multiplication function to calculate μ i . The MODEL statement specifies the likelihood function as given in Equation 2. The first step in evaluating the results is to review the convergence diagnostics. With ODS Graphics turned on, PROC MCMC produces graphs. Figure 2 displays convergence diagnostic graphs for the β 0  regression parameter. The trace plot indicates that the chain does not appear to have reached a stationary distribution and appears to have poor mixing. The diagnostic plots for the rest of the parameters (not shown here) tell a similar story.   Figure 2 Bayesian Diagnostic Plots for β 0     The non-convergence exhibited here results because the parameters are scaled very differently from each other for these data. The random walk Metropolis algorithm is not an optimal sampling algorithm in the case where the parameters have vastly different scales. Standardized covariates (Mayer and Younger; 1976) eliminate this problem, and the random walk Metropolis algorithm proceeds smoothly.   Bayesian Linear Regression Model with Standardized Covariates   Suppose you want to fit the same Bayesian linear regression model, but you want to use standardized covariates. You rewrite the mean function in Equation 1 as where X* is the design matrix constructed from a column of 1s and p standardized covariates. The regression parameters on the standardized scale are represented by β*. The standardized covariates are computed as follows:   for i = 1 ,..., n players and j = 1,..., p covariates, and where m j   and s j  are the mean and standard deviation of the j th covariate, respectively. The following statements manipulate the SUM_BASEBALL output data set from the earlier use of PROC MEANS. The statements create macro variables for the means and standard deviations to use later in the analysis. The macro variables are independent of SAS data set variables and can be referenced in SAS procedures to facilitate computations. The TRANSPOSE procedure transposes the SUM_BASEBALL data set and a DATA step creates the macro variables by using the SYMPUTX functions. The %PUT statements enable you to verify that the macro variables have been created successfully. proc transpose data=sum_baseball out=tab; id _stat_; run; data _null_; set tab; sub = put((_n_-1), 1.); call symputx(compress('m' || sub,'*'), mean); call symputx(compress('s' || sub,'*'), std); run; %put &m1 &m2 &m3 &m4 &m5 &m6 &m7 &m8; %put &s1 &s2 &s3 &s4 &s5 &s6 &s7 &s8;   In this example, m j   and s j  were calculated in the MEANS procedure and recorded in the macro variables M1–M8 and S1–S8, respectively. The STANDARD procedure computes standardized values of the variables in the original data set. proc standard data=baseball out=baseball_std mean=0 std=1; var no_hits -- cr_hits2; run;   The new likelihood function for the logarithm of the salary and corresponding standardized covariates is as follows:   =    For ease of interpretation and inference, you can transform the standardized regression parameters back to the original scale with the following formulas:   Suppose the following diffuse prior distribution is placed on β* : The prior distribution for σ 2  is given in Equation 3. Using Bayes' theorem, the likelihood function and prior distributions determine the posterior distribution of the parameters as follows:     The following SAS statements fit the Bayesian linear regression model. The MONITOR= option outputs analysis on selected symbols of interest in the program. ods graphics on; proc mcmc data=baseball_std nbi=10000 nmc=20000 seed=1181 propcov=quanew monitor=(beta0-beta8 sig2); array beta[9] beta0-beta8 (0); array betastar[9] betastar0-betastar8; array data[9] 1 no_hits no_runs no_rbi no_bb yr_major cr_hits yr_major2 cr_hits2; array mn[9] (0 &m1 &m2 &m3 &m4 &m5 &m6 &m7 &m8); array std[9] (0 &s1 &s2 &s3 &s4 &s5 &s6 &s7 &s8); parms betastar: 0; parms sig2 1; prior betastar: ~ normal(0,var = 1000); prior sig2 ~ igamma(shape = 3/10, scale = 10/3); call mult(betastar, data, mu); model logsalary ~ n(mu, var = sig2); beginprior; summ = 0; do i = 2 to 9; beta[i] = betastar[i]/std[i]; summ = summ + beta[i]*mn[i]; end; beta0 = betastar0 - summ; endprior; run; ods graphics off;   The first two ARRAY statements specify names for the regression coefficients β and β* for the original and standardized scale, respectively. The last three ARRAY statements for DATA, MN, and STD vectors take advantage of PROC MCMC’s ability to use both matrix functions and the SAS programming language. The PARMS, PRIOR, and MODEL statements are called with the same syntax as in the first call to the MCMC procedure. The BEGINPRIOR and ENDPRIOR statements reduce unnecessary observation-level computations. The statements inside the BEGINPRIOR and ENDPRIOR statements create a block of statements that are run only once per iteration rather than once for each observation at each iteration. This enables a quick update of the symbols enclosed in the statements. The statements within the BEGINPRIOR and ENDPRIOR block transform the β* sampled values back to β. The trace plot in Figure 3 indicates that the chain appears to have reached a stationary distribution. It also has good mixing and is dense. The autocorrelation plot indicates low autocorrelation and efficient sampling. Finally, the kernel density plot shows the smooth, unimodal shape of posterior marginal distribution for β 0 . The remaining diagnostic plots (not shown here) similarly indicate good convergence in the other parameters. Using standardized covariates solves the case of convergence for this model for these data.   Figure 3 Bayesian Diagnostic Plots for β 0       Figure 4 reports summary and interval statistics of all parameters. For example, the mean salary increases by an estimated factor of   (approximately 27%) for each year the player was in Major League Baseball. Similarly, using the same formula,  , you can see how the mean salary changes by one unit for each of the p covariates. Both the equal tail and highest posterior density (HPD) intervals include 0 forβ 1 , β 2 , and β 3 , indicating that the change in salary with respect to these covariates is not significant. The number of years played seems to be the most influential covariate, followed by the number of career hits.   Figure 4 Posterior Model Summary of Bayesian Linear Regression with Standardized Covariates   The MCMC Procedure   Posterior Summaries Parameter N Mean Standard Deviation Percentiles 25% 50% 75% beta0 20000 1.6465 0.0666 1.6006 1.6470 1.6935 beta1 20000 -0.00007 0.000938 -0.00071 -0.00004 0.000594 beta2 20000 0.000882 0.00167 -0.00023 0.000860 0.00200 beta3 20000 0.00186 0.000993 0.00119 0.00186 0.00253 beta4 20000 0.00218 0.000980 0.00152 0.00217 0.00281 beta5 20000 0.1042 0.0205 0.0902 0.1038 0.1176 beta6 20000 0.000748 0.000163 0.000642 0.000750 0.000857 beta7 20000 -0.00629 0.000978 -0.00692 -0.00629 -0.00562 beta8 20000 -1.46E-7 5.867E-8 -1.86E-7 -1.47E-7 -1.08E-7 sig2 20000 0.0595 0.00533 0.0558 0.0592 0.0629     Posterior Intervals Parameter Alpha Equal-Tail Interval HPD Interval beta0 0.050 1.5123 1.7714 1.5120 1.7705 beta1 0.050 -0.00195 0.00165 -0.00192 0.00168 beta2 0.050 -0.00235 0.00417 -0.00233 0.00418 beta3 0.050 -0.00006 0.00382 -0.00001 0.00383 beta4 0.050 0.000236 0.00412 0.000303 0.00416 beta5 0.050 0.0651 0.1450 0.0625 0.1415 beta6 0.050 0.000428 0.00107 0.000428 0.00107 beta7 0.050 -0.00827 -0.00443 -0.00822 -0.00442 beta8 0.050 -2.62E-7 -3.17E-8 -2.63E-7 -3.32E-8 sig2 0.050 0.0498 0.0705 0.0494 0.0699   References   Mayer, L. S. and Younger, M. S. (1976), “Estimation of Standardized Regression Coefficients,” Journal of the American Statistical Association , 71(353), 154–157.   Reichler, J. L., ed. (1987), The 1987 Baseball Encyclopedia Update , New York: Macmillan.   Time Inc. (1987), “What They Make,” Sports Illustrated , 54–81. ... View more

Overview   Health service researchers frequently study length of hospital stay (LOS) as a health outcome. Generally originating from heavily skewed distributions, LOS data can be difficult to model with a single parametric model. Mixture models can be quite effective in dealing with such data. This example illustrates how to perform a Bayesian analysis of an exponential mixture model for LOS data. The experimental MCMC procedure is used for this analysis. Analysis   The LOS data analyzed in this example originate from geriatric patients in a psychiatric hospital in North East London in 1991 and were studied by Harrison and Millard ( 1991 ) and McClean and Millard ( 1993 ) . Each observation represents the LOS in days for an admitted patient. data inputdata; input los @@; datalines; 1671 1300 722 586 552 525 364 359 321 302 272 248 226 216 208 182 141 141 132 120 117 115 114 113 104 103 101 99 96 94 93 92 88 84 83 81 79 74 70 63 62 62 61 56 55 53 53 51 51 50 49 36 33 33 33 29 28 26 24 19 16 16 15 ... more lines ... 35 28 20 19 9 5 2 2 2 22317 14006 11549 11006 8981 8402 7947 7266 6693 4408 4010 4003 1970 1857 1849 1833 1770 1769 1514 1217 956 924 611 386 280 93 ; Before using this data in a mixture model setting, use the MEANS procedure to obtain summary statistics for the LOS data in table form. Figure 1 displays summary statistics for this data.   Figure 1 PROC MEANS Summary of LOS   The MEANS Procedure   Analysis Variable : los N Mean Std Dev Minimum Maximum 469 3712.36 5675.97 1.0000000 24028.00   Figure 2 displays a plot that illustrates the skewness of the data. The histogram of the data is overlaid with a scaled exponential density, and you can see that a single exponential density does not fit the lower values of LOS well.   Figure 2 Length of Stay, Exponential Density     Bayesian Mixture of Exponential Model   Mixture models arise naturally when one mechanism generates data according to one model and another mechanism generates data according to a different model. In this data set, the variable that indicates which mechanism generated an observation is not recorded, and only the response variable is available. A mixture model is useful for analyzing this data set because it unveils the latent heterogeneity that arises from a latent categorical variable (Fruhwirth-Schnatter; 2006 ) . Mixture models can be expressed as a weighted average of K component densities. More formally, Y is said to arise from a finite mixture model with probability density function   where  For all  ,   is the component probability density function, and   are the weights defined by the following constraints:    and  . Congdon ( 2003 ) states that exponential mixture models with relatively small number of components are effective in modeling skewed LOS data. You can write a two-component exponential mixture model for LOS data with density as follows:   for patients i = 1,..., 469. There are three parameters in the density: B, D, and 𝜂. Researchers Harrison and Millard ( 1991 ) suggest the mean parameters,B and D, as the average LOS for the standard- and long-stay groups, respectively. In addition, 𝜂 represents an unknown fraction of patients in the standard-stay group with  . Suppose the following prior distributions are placed on the three parameters:     where    indicates a prior distribution and   is the density function for the inverse-gamma distribution. Priors of this type are often called diffuse priors . The uniform (0, 1) prior expresses your lack of knowledge about the mixture proportion. Label-switching is a common problem that arises in mixture models. It was described by as a result of the invariance of the mixture likelihood function under the relabeling of the mixture components. In an effort to remove label-switching, Harrison and Millard ( 1991 ) place an identifiability constraint that the average LOS of the standard-stay group is less than that of the long-stay group; that is, B < D. Using Bayes’ theorem, the likelihood function and prior distributions determine the posterior distribution of B, D, and π as follows:   PROC MCMC obtains samples from the desired posterior distribution, which is determined by the prior and likelihood specified. It does not require the form of the posterior distribution. The following SAS statements use the prior distributions to fit the Bayesian exponential mixture model. The PROC MCMC statement invokes the procedure and specifies the input data set. The NMC= option specifies the number of posterior simulation iterations. The THIN=5 option specifies that one of every five samples is kept. The PROPCOV=QUANEW option uses the estimated inverse Hessian matrix as the initial proposal covariance matrix. ods graphics on; proc mcmc data=inputdata seed=1010 nmc=50000 thin=5 propcov=quanew; ods output PostSummaries = post_summ; parms B 100 D 6000 pi 0.5; prior B D ~ igamma(3/10, scale = 10/3); prior pi ~ uniform(0,1); if (B < D) then llike = log(pi*pdf("expo", los, B) + (1-pi)*pdf("expo", los, D)); else llike = .; model general(llike); run; ods graphics off; The ODS OUTPUT statement creates an output data set post_summ used later for a graphical analysis of the fit in Figure 8. The PARMS statement puts all three parameters B, D, and 𝜂. in a single block and assigns initial values to each of them. The PRIOR statements specify priors for all the parameters. Note that B and D can be specified with one PRIOR statement because they have the same prior distribution. The IF-ELSE statements enable different values of LOS to have different log-likelihood functions, depending on whether the order constraint placed on the mean parameters is satisfied. The MODEL statement specifies that llike is the log likelihood for each observation in the model and is simply a missing value when the order constraint is not met. By turning ODS Graphics on, PROC MCMC produces graphs at the end of the procedure which enable you to visually examine the convergence of the chain. See Figure 3. Inferences should not be made if the Markov chain has not converged.   Figure 3 LOS Diagnostic Plots for B, D, and 𝜂.      Figure 3 displays convergence diagnostic graphs for parameters B, D, and 𝜂.The trace plots indicate that the chains appear to have reached stationary distributions. The chain also has good mixing and is dense. The autocorrelation plots indicate low autocorrelation and efficient sampling. Finally, the kernel density plots show the smooth, unimodal shape of posterior marginal distributions for each parameter. If label-switching had occured, you would see jumps in the trace plots or multimodal density plots. Figure 4 contains the "Parameters" table which lists the names of the parameters, the blocking information, the sampling method used, the starting values, and the prior distributions.   Figure 4 LOS Parameter Information   Parameters Parameter Sampling Method Initial Value Prior Distribution B N-Metropolis 100.0 igamma(3/10, scale = 10/3) D N-Metropolis 6000.0 igamma(3/10, scale = 10/3) pi N-Metropolis 0.5000 uniform(0,1)   The "Tuning History" table, shown in Figure 5 , displays how the tuning stage progresses for the multivariate random walk Metropolis algorithm used by PROC MCMC to generate samples from the posterior distribution. An important aspect of the algorithm is the calibration of the proposal distribution. The tuning of the Markov chain is broken into a number of phases. In each phase, PROC MCMC generates trial samples and automatically modifies the proposal distribution as a result of the acceptance rate. The "Burn-In History" and the "Sampling History" tables show the burn-in and main phase sampling, respectively.   Figure 5 LOS Burn-In and Sampling History   Tuning History Phase Block Scale Acceptance Rate 1 1 2.3800 0.0280 2 1 1.0123 0.0700 3 1 0.5221 0.3420   Burn-In History Block Scale Acceptance Rate 1 0.5221 0.3980   Sampling History Block Scale Acceptance Rate 1 0.5221 0.3811   Figure 6 displays summary and interval statistics for each posterior distribution.   Figure 6 LOS MCMC Summary and Interval Statistics The MCMC Procedure   Posterior Summaries Parameter N Mean Standard Deviation Percentiles 25% 50% 75% B 10000 619.1 74.2261 568.9 618.3 668.3 D 10000 7752.6 721.2 7267.1 7710.4 8202.1 pi 10000 0.5638 0.0403 0.5370 0.5655 0.5916   Posterior Intervals Parameter Alpha Equal-Tail Interval HPD Interval B 0.050 473.6 767.0 472.0 763.7 D 0.050 6418.1 9288.4 6339.5 9195.0 pi 0.050 0.4826 0.6389 0.4825 0.6387   Figure 7 reports a number of convergence diagnostics to assist in determining convergence. These are the Monte Carlo standard errors, the autocorrelations at selected lags, the Geweke diagnostics, and the effective sample sizes.   Figure 7 LOS MCMC Convergence Diagnostics The MCMC Procedure   Monte Carlo Standard Errors Parameter MCSE Standard Deviation MCSE/SD B 3.9598 74.2261 0.0533 D 38.2006 721.2 0.0530 pi 0.00164 0.0403 0.0408   Posterior Autocorrelations Parameter Lag 1 Lag 5 Lag 10 Lag 50 B 0.9150 0.6588 0.4459 0.0689 D 0.9220 0.6766 0.4760 0.0354 pi 0.5086 0.3438 0.2564 0.0267   Geweke Diagnostics Parameter z Pr > |z| B 0.1718 0.8636 D -0.2669 0.7895 pi -0.3073 0.7586   Effective Sample Sizes Parameter ESS Correlation Time Efficiency B 351.4 28.4595 0.0351 D 356.4 28.0580 0.0356 pi 601.9 16.6129 0.0602   Figure 6  displays the marginal posterior summaries. The larger group of standard-stay patients has an average LOS of B = 619 days relative to the smaller group of long-stay patients whose stay averages D = 7752 days. The posterior average mixture proportions are 56% and 44% for the standard-stay and long-stay group, respectively. Figure 8 illustrates the scaled mixture of exponential density and the gain in model fit from the two-component exponential mixture model.   Figure 8 Length of Stay, Exponential and Mixture Density     The single-component exponential model had an average LOS of 3712 days. The mean parameters found when fitting an exponential mixture model to the standard-stay group and the long-stay group are 619 and 7752 days, respectively. The mixture distribution fits the data better than the exponential distribution, especially at the low values of LOS. Additional components could be fit and evaluated in similar methods or with information criteria.   References   Congdon, P. (2003), Applied Bayesian Modeling , John Wiley & Sons.   Fruhwirth-Schnatter, S. (2006), Finite Mixture and Markov Switching Models , Springer. Harrison, G. and Millard, P. (1991), “Balancing Acute and Long-Term Care: The Mathematics of Throughput in Departments of Geriatric Medicine,” Meth. Infor. Medicine , 30, 221–228.   McClean, S. and Millard, P. (1993), “Patterns of Length of Stay after Admission in Geriatric Medicine: An Event History Approach,” The Statistician , 42(3), 263–274. ... View more

Overview   This example illustrates fitting Bayesian zero-inflated Poisson (ZIP) models to zero-inflated count data with the experimental MCMC procedure. ZIP models are often used when count data show an excess number of zeros, which in turn causes over-dispersion. Consider survey data collected at a state park concerning the number of fish that visitors had caught in the last six months. People from two different populations provided information: those who attended the park and fished, and those who attended the park and did not fish. Zero fish caught has a different meaning for the two populations.   Analysis   Count data frequently display over-dispersion and excess zeros, which motivates zero-inflated count models (Lambert; 1992; Greene; 1994) . Zero-inflated count models offer a way of modeling the excess zeros in addition to allowing for over-dispersion in a standard parametric model. Zero inflation arises when one mechanism generates only zeros and the other process generates both zero and nonzero counts. Zero-inflated models can be expressed as a two-component mixture model where one component has a degenerate distribution at zero and the other is a count model. More formally, a zero-inflated model can be written as   =   where      (1)     and where     is a standard count model with mean  , support  , and 𝜂 is a mixture proportion with  .   The following hypothetical data represents the number of fish caught by visitors at a state park. Variables are created for the visitor’s age and gender. Two dummy variables, FEMALE and MALE, are created to indicate the gender. data catch; input gender $ age count @@; if gender = 'F' then do; female = 1; male = 0; end; else do; female = 0; male = 1; end; obs = _N_; datalines; F 54 18 M 37 0 F 48 12 M 27 0 M 55 0 M 32 0 F 49 12 F 45 11 M 39 0 F 34 1 F 50 0 M 52 4 M 33 0 M 32 0 F 23 1 F 17 0 F 44 5 M 44 0 F 26 0 F 30 0 F 38 0 F 38 0 F 52 18 M 23 1 F 23 0 M 32 0 F 33 3 M 26 0 F 46 8 M 45 5 M 51 10 F 48 5 F 31 2 F 25 1 M 22 0 M 41 0 M 19 0 M 23 0 M 31 1 M 17 0 F 21 0 F 44 7 M 28 0 M 47 3 M 23 0 F 29 3 F 24 0 M 34 1 F 19 0 F 35 2 M 39 0 M 43 6 ; Although these data appear to be a likely candidate for a ZIP model, you typically begin with a standard analysis and evaluate the evidence for over-dispersion.   Bayesian Poisson Regression Model   Suppose you want to fit a Bayesian Poisson regression model for the number of fish caught with density as follows:   (2)      for the i = 1,....,52 surveyed park visitors. The likelihood function for each of the counts and corresponding covariates is   (3)   where   denotes a conditional probability density. The Poisson density is evaluated at the specified value of   and corresponding mean parameter  . The three parameters in the likelihood are  , and  , which correspond to an intercept, slope for age for females, and slope for age for males, respectively. Suppose the following prior distributions are placed on the three parameters:     (4)   where   indicates a prior distribution. The diffuse   prior expresses your lack of knowledge about the regression parameters. Using Bayes’ theorem, the likelihood function and prior distributions determine the posterior distribution of  , and   as follows:   PROC MCMC obtains samples from the desired posterior distribution, which is determined by the prior and likelihood specified. It does not require the form of the posterior distribution. The goodness-of-fit Pearson chi-square statistic   as given in McCullagh and Nelder (1989) , is calculated to assess model fit:   (5)     for a Poisson likelihood where    with   as defined in Equation 2,    is the expectation, and   is the variance. If there is no over-dispersion, the Pearson statistic would roughly equal the number of observations in the data set minus the number of parameters in the model. The following SAS statements use the diffuse prior distributions to fit the Bayesian Poisson regression model and calculate the fit statistic. ods graphics on; proc mcmc data=catch seed=1181 nmc=100000 thin=10 propcov=quanew monitor =(_parms_ Pearson); ods select Parameters PostSummaries PostIntervals tadpanel; parms beta0 0 beta1 0 beta2 0; prior beta: ~ normal(0,var=1000); mu = exp(beta0 + beta1*female*age + beta2*male*age); model count ~ poisson(mu); if obs = 1 then Pearson = 0; Pearson = Pearson + ((count - mu)**2/mu); run; ods graphics off;   The PROC MCMC statement invokes the procedure and specifies the input data set. The SEED= option specifies a seed for the random number generator, which guarantees the reproducibility of the random stream. The NMC= option specifies the number of posterior simulation iterations. The THIN=10 option controls the thinning of the Markov chain and specifies that one of every 10 samples is kept. The PROPCOV=QUANEW option initializes the Markov chain at the posterior mode and uses the estimated inverse Hessian matrix as the initial proposal covariance matrix in the random walk Metropolis algorithm. The MONITOR= option specifies a list of symbols (which can be either parameters or functions of the parameters in the model) for which inference is to be done. The symbol _parms_ is a shorthand for all model parameters—in this case,  , and  . The symbol Pearson refers to the Pearson chi-square statistic, which is a function of the data and parameters. The PARMS statement puts all three parameters  , and   in a single block and assigns initial values to each of them. The PRIOR statement specifies priors for all the parameters. The notation beta: in the PRIOR statement is a shorthand for all variables that start with ' beta '. In this example, it includes beta0 , beta1 , and beta2 . The shorthand notation is not necessary, but it keeps your code succinct. The assignment statement for mu calculates   in the Poisson model, as given in Equation 2. The MODEL statement specifies the likelihood function for COUNT, as given in Equation 3. The next two lines of statements use the posterior samples of mu and the data set variable count to calculate the Pearson chi-square statistic. The IF statement resets the value of Pearson to be zero at the top of the data set (that is, when the data set variable obs is 1). As PROC MCMC cycles through the data set at each iteration, the procedure cumulatively adds the Pearson chi-square statistic over each value of count . By the end of the data set, you obtain the Pearson chi-square statistic, as defined in Equation 5. It is essential to examine the convergence of the Markov chains before you proceed with posterior inference in Bayesian analysis. With ODS Graphics turned on, PROC MCMC produces graphs at the end of the procedure output which allow you to visually examine the convergence of the chain. See Figure 1. Inferences should not be made if the Markov chain has not converged.   Figure 1 Bayesian Poisson Model Diagnostic Plots for  , and      Figure 1 displays convergence diagnostic plots for  , and  . The trace plots show that the mean of the Markov chain is constant over the graph and is stabilized. The chain was able to traverse the support of the target distribution, and the mixing is good. The trace plots show that the Markov chain appears to have reached stationary distributions. The autocorrelation plots indicate low autocorrelation and efficient sampling. The kernel density plots show smooth, unimodal posterior marginal distributions for each parameter. PROC MCMC produces formal diagnostic tests by default, but they are omitted here because an informal check on the chains, autocorrelation, and posterior density plots show desired stabilization and convergence. The "Parameters" tables, shown in Figure 2, lists the names of the parameters, the sampling method used, the starting values, and the prior distributions.   Figure 2 Bayesian Model Information The MCMC Procedure   Parameters Parameter Sampling Method Initial Value Prior Distribution beta0 N-Metropolis 0 normal(0,var=1000) beta1 N-Metropolis 0 normal(0,var=1000) beta2 N-Metropolis 0 normal(0,var=1000)   Figure 3 displays summary and interval statistics for each parameter’s posterior distribution. PROC MCMC also calculates the sampled value of the Pearson chi-square at each iteration and produces posterior summary statistics for it.   Figure 3 Posterior Model Summary of Poisson Regression   The MCMC Procedure   Posterior Summaries Parameter N Mean Standard Deviation Percentiles 25% 50% 75% beta0 10000 -4.0186 0.5551 -4.3787 -4.0024 -3.6411 beta1 10000 0.1284 0.0117 0.1204 0.1281 0.1361 beta2 10000 0.1047 0.0125 0.0963 0.1043 0.1129 Pearson 10000 92.2674 12.2613 83.3871 90.0545 98.8112   Posterior Intervals Parameter Alpha Equal-Tail Interval HPD Interval beta0 0.050 -5.1475 -2.9799 -5.0903 -2.9400 beta1 0.050 0.1064 0.1518 0.1065 0.1518 beta2 0.050 0.0808 0.1295 0.0799 0.1283 Pearson 0.050 74.8858 121.7 72.8727 116.1   With n = 52 and three model parameters, the sampled value 92.2674 of the Pearson chi-square statistic is much greater than 52 - 3 = 49   , providing evidence of over-dispersion. Zero inflation is a likely cause of this over-dispersion. A Bayesian ZIP model accounts for the extra zeros and potentially provides a better fit to the data.   Bayesian ZIP Regression Model   You can write a Bayesian ZIP regression model for the number of fish caught as follows:     where   is defined in Equation 1,   is defined in Equation 2, and   for the i = 52 surveyed visitors. The model is a weighted average of the degenerate function (which places all mass at zero) and the Poisson regression.   The likelihood function for each of the counts and corresponding covariates is      (6)   where   denotes a conditional probability density and the Poisson density is evaluated at the specified value of   and corresponding mean parameter  . The degenerate distribution   is one when COUNT equals zero and remains zero for a COUNT value greater than zero. The four parameters in the likelihood are  , and  , which correspond to an intercept, slope for age for females, slope for age for males, and the mixture proportion, respectively. Suppose again that the three regression parameters have the same diffuse, normal priors as in Equation 4. The mixture proportion has a uniform(0,1) prior distribution. The Pearson chi-square statistic for the ZIP model is calculated as in Equation 5, but now the mean and the variance respectively are   The following SAS statements use the prior distributions to fit the Bayesian ZIP regression model and calculate the Pearson chi-square statistic. ods graphics on; proc mcmc data=catch seed =1181 nmc=100000 thin=10 propcov=quanew monitor =(_parms_ Pearson); ods select Parameters PostSummaries PostIntervals tadpanel; parms beta0 0 beta1 0 beta2 0 eta .3; prior beta: ~ normal(0,var=1000); prior eta ~ uniform(0,1); mu=exp(beta0 + beta1*female*age + beta2*male*age); llike=log(eta*(count eq 0) + (1-eta)*pdf("poisson",count,mu)); model general(llike); if obs = 1 then Pearson = 0; mean = (1 - eta)*mu; sigma2 = (1 - eta)*mu*(1 + eta*mu); Pearson = Pearson + ((count - mean)**2/sigma2); run; ods graphics off;   The parameter 𝜂 and its starting value are added to the PARMS statement with the regression parameters and their starting values. The PRIOR statement remains the same for the regression parameters, but an additional PRIOR statement is needed for the mixture proportion. The assignment statement for mu calculates the expected value of COUNT in the Poisson model, as given in Equation 2. The assignment statement for llike evaluates the log density of Equation 6. The expression (count eq 0) in llike acts as an indicator variable for the degenerate distribution   ; it is one when COUNT equals zero and zero for values of COUNT greater than zero. The MODEL statement specifies that llike is the log likelihood for each observation in the model. The Pearson chi-square statistic is calculated according to Equation 5; the moments are evaluated in the mean and sigma2 assignment variables.   Figure 4 Bayesian ZIP Diagnostic Plots for  ,  , and 𝜂            The diagnostic plots for the regression parameters and mixture proportion are illustrated in Figure 4. They show the desired convergence, low autocorrelation, and smooth unimodal marginal posterior densities for the parameters. Figure 5 displays the "Parameters" table for the ZIP regression model. The "Parameters" table now includes information about the mixture proportion parameter.   Figure 5 Bayesian ZIP Regression Model Information   The MCMC Procedure   Parameters Parameter Sampling Method Initial Value Prior Distribution beta0 N-Metropolis 0 normal(0,var=1000) beta1 N-Metropolis 0 normal(0,var=1000) beta2 N-Metropolis 0 normal(0,var=1000) eta N-Metropolis 0.3000 uniform(0,1)     Figure 6 displays summary and interval statistics for each parameter’s posterior distribution. PROC MCMC also calculates the sampled value of the Pearson chi-square at each iteration. This statistic in the Bayesian ZIP model is greatly reduced with a value of 48.8650, which suggest a much better fit compared to its value of 92.2674 in the Poisson model.   Figure 6 Posterior Summary of Bayesian ZIP Regression   The MCMC Procedure   Posterior Summaries Parameter N Mean Standard Deviation Percentiles 25% 50% 75% beta0 10000 -3.5345 0.6349 -3.9516 -3.5278 -3.0980 beta1 10000 0.1217 0.0132 0.1127 0.1214 0.1304 beta2 10000 0.1054 0.0138 0.0961 0.1051 0.1146 eta 10000 0.3074 0.0936 0.2420 0.3050 0.3719 Pearson 10000 48.8650 10.2696 41.4549 47.2985 54.4092   Posterior Intervals Parameter Alpha Equal-Tail Interval HPD Interval beta0 0.050 -4.8193 -2.3170 -4.7908 -2.2925 beta1 0.050 0.0963 0.1485 0.0955 0.1476 beta2 0.050 0.0786 0.1331 0.0785 0.1328 eta 0.050 0.1296 0.4926 0.1290 0.4915 Pearson 0.050 33.7757 73.5911 31.6653 68.8728   Modeling the catch data set with a Bayesian ZIP regression model accounts for the zero inflation and removes the overdispersion in the Poisson regression model. Figure 6 shows the posterior parameter summaries in addition to the lowered Pearson chi-square statistic. The posterior mean of the mixture probability is 0.3074 for the zero-inflated component and 1 - 0.3074 = 0.6926 for the Poisson regression component. The posterior parameter means of β 1 and β 2   are 0.1217 and 0.1054, respectively. That is, an increase in one year of age is estimated to be associated with a change in the mean number of fish caught by females and males by a factor of exp(β 1 ) = 1.13 and exp(β 2 ) = 1.1, respectively. For this model, choice of priors, and data set, females catch more fish with age than males.   References   Greene, W. H. (1994), Accounting for Excess Zeros and Sample Selection in Poisson and Negative Binomial Regression Models , working paper 94-10, New York University, Leonard N. Stern School of Business, Department of Economics, available at http://ideas.repec.org/p/ste/nystbu/94-10.html .   Lambert, D. (1992), “Zero-Inflated Poisson Regression Models with an Application to Defects in Manufacturing,” Technometrics , 34, 1–14.   McCullagh, P. and Nelder, J. A. (1989), Generalized Linear Models , Second Edition, London: Chapman & Hall. ... View more

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