Evidence for Autoregressive Conditional Heteroskedastic errors in growth traits of beef cattle

Evidence is presented for ‘generalized autoregressive conditional heteroskedasticity’ processes (GARCH(1,1)), in the residuals of beef cattle growth traits. This process can account for differences in variance at different time points, with the advantage of using a parsimonious parametrization. Data used were 10271 birth weights (BW), 19992 weaning weights (WW) and 9717 weight at 18 months (FW), from five herds registered in the national evaluation of the Brangus breed in Argentina. The residuals calculated from the 2005 genetic evaluation were regressed on Julian dates by least squares. From a second set of residuals out of the linear regression model, Maximum Likelihood estimation via the Fisher scoring algorithm was used to estimate the GARCH(1,1) parameters. Eight out of fifteen one-sided Lagrange multiplier statistics significantly (P < 0.05) rejected the hypothesis of null GARCH(1,1) parameters in the genetic evaluation residuals. Incorporating these effects in genetic evaluation is feasible due to the diagonal covariance matrix induced by the process on each trait, which simplifies building the mixed model equations.


Introduction
The presence of heterogeneous error variances in field data used for genetic evaluation has been extensively documented in different species and different traits.Fail to account for heterogeneous variance in models of genetic evaluation reduces intensity and accuracy of selection (HILL, 1984).The reason is that animals displaying higher performances in highly variable herds tend to be selected over individuals that are high performers in herds with low variability (HILL, 1984; VAN VLECK, 1987).As a consequence, genetic progress will be reduced.On the other hand, if estimates of the heterogeneous variances are available, the mixed model equations for Best Linear Unbiased Prediction (HENDERSON, 1984) can be modified to account for heteroskedasticity into the model (GIANOLA, 1986).However, several covariance structures that account for heterogeneity of variances and covariances are either not amenable to computation (i.e.non-diagonal R covariance matrices of error effects), or are difficult to estimate due to the high number of parameters involved.J. L. Foulley andcoworkers (FOULLEY et al, 1990, 1992;SAN CRISTOBAL et al., 1993) have developed models to detect causes of heteroskedasticity based on taking natural logs of residual variances and setting a linear model for the transformed dispersion parameters.There may be instances in which the existing heterogeneity is time-dependent, and can not be ascribed to a single cause, such as contemporary group, age of dam, sex, etc.In these cases, a stochastic process may be envisaged at the herd level that makes the variance of observations from animals that are born and raised one or more days apart, as related to each other.The process is such that its unconditional variance does not display the heterogeneity, but the variance conditional on one or more previous errors does.Within the econometrics literature, ENGLE (1982) has described such a process, and called it 'autoregressive conditional heteroskedasticty'.The literature is abundant in related processes such as GARCH, IGARCH, ARCH-M, to quote some.ARCH-GARCH processes can take into account periods of low variance and periods where the variability observed is high.The objective of this research is to display evidence for the presence of GARCH processes in growth traits of beef cattle at the herd level.

Data
Growth data used in the study consisted of 10271 birth weights (BW), 19992 weaning weights (WW) and 9717 weights at approximately 18 months of age (FW), of bulls and heifers registered in the animal genetic evaluation program (ERBra) of the Argentine Brangus Association.Out of the 56 herds that participated of ERBra 2005, five were chosen for the analysis.The criteria used in selecting the herds were 1) number of years of recording, 2) completeness of the data, and 3) regional representativity.The day in which the first recorded animal was born was arbitrarily set to day 1 for each herd, and the day of birth (t) of any animal born later was scored in Julian days.Residuals for BW, WW and FW, were calculated from the multiple trait animal model evaluation used in the ERBra.Fixed effects included contemporary groups for all traits, sex of calf and age of dam for BW and WW, and age of weaning for WW and age at the day of measure for FW.In the ERBra, the parametrization of HILL (1982) is used to correct for heterosis in WW using mean dominance and mean additive maternal effects, whereas dominance is included in the model for FW.Random direct breeding values were included for the three traits, whereas WW also included maternal breeding values.Estimates of all fixed effects plus predicted breeding values from the 2005 evaluation were employed to obtain the estimated residuals.

Estimation and hypothesis testing
The characteristics and properties of ARCH and GARCH processes are described in Appendix A. The parameters of the GARCH(1,1) processes were estimated by maximum likelihood using the Fisher scoring algorithm.The estimating equations for the ARCH(1) process are described in Appendix B, whereas the derivatives needed for the GARCH(1,1) process are obtained in Appendix C. In order to test for the presence of a GARCH (1,1) processes in the data, we used the Likelihood Multiplier (LM) statistics originally described by ENGLE (1982) to test for an ARCH(1) process.LEE (1991) showed that the LM test for an ARCH(1) is numerically identical to the LM test for the GARCH(1,1).Therefore, we used the LM statistics to test directly for the GARCH(1,1), where the null hypothesis is " 1 = 0 and $ 1 = 0.A further refinement consisted in using the more powerful one-sided version of the LM test discussed by DEMOS and SENTANA (1998), which will be explained below.
In the context of beef cattle growth traits, the LM procedure consists on fitting a linear regression model on the residuals from the genetic evaluation on julian dates, and then test for the presence of a GARCH process in the residuals from this regression.A problem with animal breeding data is that calvings occur irregularly through time, so that in any given date there may be several animals being born.As standard econometric techniques are based on fixed time interval analysis, ENGLE and RUSELL (1998) expanded ARCH-GARCH processes to account for variable time interval analysis into autoregressive conditional duration (ACD) models.However, ENGLE (2002) observed that standard ARCH-GARCH techniques can be used after transforming the ACD irregularly spaced random variable x t (in our case Day of Birth, DOB, in Julian days) to t x , and then setting the mean to zero.Let e (ij)t be the residual of trait j from animal i with transformed DOB (denoted as t).The residual e (ij)t was regressed on the transformed DOB of individual i (t i ) such that: where * 0 is the intercept and * 1 is the slope of the residuals through time, and ( ) ε is the error term.It is expected that the 'regression residuals' from the genetic evaluation do not deviate from a zero expectation, i.e. * 0 = 0 and * 1 = 0.Moreover, if there is no evidence for an ARCH or GARCH process, then the term ( ) ε is just 'white noise': normal independent random variables with zero mean and constant variance.Regression residuals were obtained with least squares estimates of the parameters in (8) as follows The LM test for a GARCH(1,1) process was implemented for each of the three traits and each of the five herds.The null hypothesis was H 0 : α 1 = 0 and $ 1 = 0, and the test statistics is equal to ( ) where F = [F 1 :..: under the null hypothesis, and . The R 2 is the squared multiple correlation between f 0 and F. Under the null hypothesis of the one-sided test proposed by DEMOS and SENTANA (1998), the test statistics is asymptotically distributed as a 50:50 mixture of c 2 (0) and c 2 (1) , which was the distribution used to calculate the p-values.

Results
Tables 1, 2, 3 display the LM tests for the hypothesis of the GARCH(1,1) process (H 0 : 1 α = 0 and ß 1 = 0) on the errors of BW, WW, and FW, respectively.There were significant (P < 0.05) evidences that 1 α > 0 and ß 1 > 0 for all three traits in herd B, and at least one trait showed a value of 1 α and ß 1 significantly (P < 0.05) greater than 0 in the other herds: WW for herd A, BW and FW for herd C, WW for herd D, and BW for herd E. Length of data collection (number of Julian dates) was not related to the result of the tests.Most values of the kurtosis for the unconditional distribution of error terms were greater than 3 whenever the LM test was significant, as expected by theory.The Figure depicts the distribution of errors for WW with time in herd A. Observe that the period from 1 to 2000 days displays larger changes in variance than the interval between 2000 and 5000 days, although the mean does not deviate from 0 throughout the 30 years of data collection.Notice the typical GARCH effect by which large (small) absolute values are followed by large (small) values of unpredictable sign.Estimates of the parameters are showed in Tables 4, 5 and 6.Across all herds and traits that displayed significant LM statistics in favor of a GARCH(1,1) process, parameter estimates ranged from 0.009 to 0.182 for α 1 , and 0.079 to 0.988 for β 1 .

Discussion
In the current research, we found evidence that the environmental variances of growth traits in beef cattle can be modeled with GARCH(1,1) processes.Heteroskedasticity in growth traits have also been reported by ROBERT-GRANIE et al (2002); RODRIGUEZ-ALMEIDA et al (1995) and SAN CRISTOBAL et al (1993).HILL (1984) observed that animals from the more variable environments tend to be selected if heterogeneous environmental variances are not accounted for.He further indicated that heterogeneity of variances induces nonlinearity between true and predicted breeding value, a fact observed on simulated data by MEUWISSEN and VAN DER WERF (1993).Thus, accuracy of prediction of breeding value will decrease due to unaccounted variance heterogeneity.There are two major advantages to the use of ARCH-GARCH processes for modeling heteroskedasticity.First, the parameters for the ARCH-GARCH process describe variance heterogeneity at the individual animal level, rather than to using groups of animals to model heteroskedasticity (say, by herd, by year, or by contemporary groups), being therefore more sensitive to differences in environmental variance.Second, the number of parameters introduced (2 per herd for ARCH(1) and 3 per herd for GARCH(1,1)) will generally be smaller than using one variance component per contemporary group.As an example, in the Argentinian Brangus evaluation the amount of information for a single parameter in the GARCH(1,1) process, is respectively, 3.3, 13.2 o 9.7 times greater for BW, WW and FW, as compared to using one variance parameter per contemporary group.We will deal with estimating the dispersion parameters in multiple-trait animal models elsewhere.If estimates are available of the ARCH or GARCH parameters, the processes can be incorporated into an animal model of genetic evaluation.Let R be the s × s environmental covariance matrix among s traits.Then, R can be decomposed as R = S C S, where S is a diagonal matrix with non-zero diagonal entries equal to the environmental standard deviations of the s traits 1 σ , …, s σ , whereas C is a a correlation matrix with diagonal elements equal to 1 and off-diagonals equal to r ij , the environmental correlation between traits i and j.To account for an ARCH-GARCH process, let R hit be the covariance matrix of environmental effects for animal i within herd h, born at time t.Now, R iht = S iht C S iht , with the diagonals of S iht equal to σ ijht : the square root of the variance in (3) for the ARCH or in (7) for the GARCH, scaled to the regular order of magnitud of R ii .
To build up the mixed model equations (MME) for a multi trait animal model, order the data by date of birth within the herd.Then, from previous solutions or a previous run, obtain the residuals for each animal within her and trait, and write them into the input data file.Next, calculate the variances σ 2 ijht while reading the data and the ijhterror term.Finally, read the data and the σ 2 ijht 's for each animal, create and invert (or obtain a g-inverse) of R ijht and form the MME as usual.
In conclusion, heteroskedasticity of growth taits in some beef herds can be modeled by GARCH(1,1) processes.Incorporating these effects in genetic evaluation is feasible due to the diagonal covariance matrix induced by the process on each trait, which simplifies building the MME.

Appendix A. The ARCH and GARCH processes
The ARCH (q) process Consider the set of random variables indexed by time T T = {e 1 ,.., e t-1 , e t ,.., e T }, and write where 0 t are standard i.i.d.normal (0,1) variables, and 2 t σ is the variance of the process at time t.The unconditional expectation of (1) is If we condition on the set T t-q , the conditional expected value is As put forward by ENGLE (1982), the variance in ( 1) is modeled conditional on T t-q with the conditions ω > 0 and 0 # α i # 1, i = 1, 2, ..., q.This equation for the variance allows capturing the dynamics of the conditional heteroskedasticity, as the conditional variance is dependent upon the realized values of e t-q ,..., e t-1 squared.In turn, this conveys the idea of the process having 'memory' on the past 1, 2.., q realizations.A further requirement is needed for the variance to be stationary, so that it does not depend on the particular time t (BOLLERSLEV et al, 1993, page 2967): ( ) Notice that the unconditional variance of (1) is constant though, and can be obtained as follows Also, the ARCH(q) variables are uncorrelated at different times: t t q t t q t t q t q t q t t q t q e e e e e e Although the e's are serially uncorrelated, they are not independent through time.In words of BOLLERSLEV et al ( 1993) "there is a tendency for large (small) absolute values of the process to be followed by other large (small) values of unpredictable sign".After (A2), ( A3) and (A4), the conditional distribution of the ARCH(q) random variable in ( 1) is equal to Note that in ARCH(q) processes the conditional probability of e 1 ,., e q-1 is unknown, as T q-1 is undefined.Although the conditional distribution of e t is time-invariant, the unconditional distribution for e t will tend to have 'fatter tails' than the normal density of , t .For example, BOLLERSLEV et al (1993) observed that the kurtosis of an ARCH(1) process with conditionally normally distributed errors is equal to: Therefore, the ARCH(1) process is leptokurtic: the curvature is high in the center of the density and the tails are 'fatter' than those of the normal distribution.
The GARCH (p, q) process BOLLERSLEV (1986) generalized ARCH(q) processes by allowing the conditional variance to be also a function on the p previous variances 2 besides the q previous realizations of e t-q ,…, e t-1 .Thus, the variance in the GARCH(p, q) process is The requirement for the variance to be stationary (HAMILTON, 1994, page 666) is as follows For the GARCH(1,1) process, BOLLERSLEV et al (1993) observed that the conditions for positivity of the conditional variance (i.e. 2 t σ > 0) are ω $ 0, α 1 $0, and Appendix B: Maximum likelihood estimation By conditioning on 1 ε from the first observation, at the start of the recursive process the likelihood does not depend on the parameters (BOLLERSLEV et al, 1993).From (A4) and (A5), the log-likelihood function can be written as By taking first derivatives of (B1) the scores are: Whereas second derivatives are equal to: The expected values of the second derivatives are needed for the Fisher scoring algorithm.Now, observe that after (1) we can write 2 2 2 t t t − = ε σ , .As the t , are N(0,1), we have that ( ) , so that the expectations of second derivatives are ( ) ( ) (1982) observed that these expectations are consistently estimated by

ENGLE
The scalars ω and ˆ1 α are the estimated dispersion parameters, and ˆt ε is the residual from the previous iteration of the Fisher scoring algorithm.The estimating equations are: successive iterations was lower than a preset value (for example, 10 -6 ) for both ω and 1 α .The equivalent system to (B6) for the GARCH(1,1) is 3 × 3 due to the added equation for 1 β .
Appendix C: Analytic derivatives of the GARCH(1,1) log-likelihood By combining expressions (3) and (4) in FIORENTINI et al (1996), the log-likelihood can be written as To take care of the first observation, Fiorentini et al (1996) employed an estimate of the unconditional expectation, so when t ≤ 0 , with the parameter vector being equal to Now, we need an expression for the derivatives of 2 t σ with respect to the elements of θ: Expression (C3) is recursive as the derivative at time t depends on to the derivatives of the variance at previous times.Thus, for t = 1, 2, 3, we have that The emerging pattern is such that In order to obtain the second derivatives, we take derivatives in (C1) to get

Figure :
Figure: Distribution of errors for weaning weight with time for herd A. (Fehlerverteilung für die Absetzgewichte abhängig von den Tagen der Herde A) t

Table 1
Results of the Lagrange multiplier tests for birth weight (Lagrangischer multipler Test für Geburtsgewichte)