Modeling Panel Time Series with Mixture Autoregressive Model

This paper considers the mixture autoregressive panel (MARP) model. This model can capture the burst and multi-modal phenomenon in some panel data sets. It also enlarges the stationarity region of the traditional AR model. An estimation method based on the EM algorithm is proposed and the assumption required of the model is quite low. To illustrate the method, we fitted the MARP model to the gray-sided voles data. Another MARP model with less restriction is also proposed.


Introduction
Panel time series data are collections of similar time series variables.A general linear dynamic model for a panel time series {X jt }, j = 1, . . ., M, t = 1, . . ., T j is given by where t denotes time and j indicate the individuals.Moreover, η t represents effects for all of the series and λ j denotes the individual impact of each of the series.
The other variables consist of the vector of explanatory variables W jt and the error term ε jt .The latter are assumed to be independent identically distributed.For details, see Hsiao (1986, p. 71) and Hjellvik and Tjφstheim (1999).Panel data have now become more and more important both in the natural science and in macroeconomic study.For instance, gray-sided voles data are yearly observations from ninety-one different positions.In macroeconomics, two very famous data sets are the National Longitudinal Survey of Labor Market Experience (NLS) and the Michigan Panel Study of Income Dynamics (PSID).Traditionally, in panel data, there are many different individuals observed through time, but the number of observations of each individual is relatively small.In many situations, time series data display non-linear features such as flat stretches, bursts and change points.At the same time, there could also be many missing datum points and outliers inside a time series data set.These are no exceptions for panel time series data.All these non-linearity and missing information make the modeling of panel time series more difficult than ordinary time series data.Many papers tried to overcome the difficulties in panel data modeling.Hjellvik and Tjφstheim (1999) treated the common factor η t in model (1.1) as nuisance parameters and estimated them by averaging the observations across all the series at fixed time point t.Fu et al. (2002b) assumed that the series were contemporaneously correlated and proposed a contemporaneous correlated threshold model for panel data analysis.In the nonparametric and semiparametric literature, Hjellvik et al. (2004) considered a nonparametric estimation and testing procedure of the nonlinearity for panel data.Lin and Ying (2001) and Fan and Li (2004) suggested to model the panel data with semiparametric method.Recently, Jin and Li (2005) considered the modeling of a contemporaneously correlated panel data with partial linear regression model.These papers partly solved the problem of change points, nonlinearity and intercorrelation.However, other features such as flat stretches, bursts and multi-modality are still untouched.Le et al. (1996) introduced the Gaussian mixture transition distribution (GMTD) model.With this model, one can capture the flat stretches, bursts and outliers in regression analysis.Wong and Li (2000) extended GMTD model to time series analysis and suggested the mixture autoregressive model (MAR) to catch the multi-modal phenomena.An N -component univariate MAR model is defined almost surely by equation (1.2).
where F (x t |F t−1 ) is the conditional cumulative distribution function of X t given the past information, evaluated at x t .The σ-field F t−1 denotes the past information set up to time t − 1 and the function Φ(•) denotes the cumulative distribution function of the standard normal distribution and the summation of the probability N i=1 π i = 1.For details, see Wong and Li (2000).From our experience, many panel time series appeared to be multi-modal.For example, Figure 1 shows four typical series in the gray-sided voles data.The four series are stationary but it is quite obvious that they are multi-modal.Hence, developing a mixture model to capture the multi-model characteristics in panel time series would be of importance.
This paper is arranged as follows.In section 2, we propose the panel mixture autoregressive model which assumes that all the panels have the same mixture probability π i , i = 1, . . ., N and one of the components has a common autoregressive structure.We develop an estimation method based on the EM algorithm.In section 3, a simulation experiment with different series length was performed to check our model and the estimation method.The minimum sample size requirement needed for the application of our model is also suggested.The gray-sided voles data are studied in section 4 as a real example of the application of our model.In section 5, we will relax the restriction that the series has the same mixture probabilities and extend the MARP model to a more flexible model.Lastly, we will draw some conclusions in section 6.

Model and Estimation
Our N -component finite mixture autoregressive model for a panel time series (MARP) X jt , j = 1, . . ., M, t = 1, . . ., T j in conditional cumulative distribution function (CDF) form, is given by where i = 1, . . ., N is the index of the components, j = 1, . . ., M is the index of the series, t = 1, . . ., T j is the time index, k is the index of the lags and p is the autoregressive order.Here, we assume that the order of each component in each series is the same.If the orders p ij are not the same, we can set p = max ij {p ij } and those φ kij = 0 when k > p ij .In the conditional CDF, the random variable X jt is evaluated at x jt given the past information up to time t − 1 which is denoted by F t−1 .The random variable X jt is taken from the ith component with probability π i where McLachlan andPeel (2000), 2000).This assumption is only for the theoretical purpose and is not used as the restriction in our model estimation.As usual, Φ(•) denotes the cumulative distribution function of the standard normal distribution.Although we assume that the order p are the same, however, the series length T j are not assumed to be same.In panel data analysis, the data often have some common features.For example, in section 4 of this paper, we studied the gray-sided voles data which were collected from different positions in Hokkaido, Japan.It is reasonable to believe that the different series would follow some common pattern.To capture the common structure among the series, we assume that the last component in each of the series has the same coefficients φ kN , k = 0, . . ., p.The error terms ε ijt are assumed to be mutually independent and from the normal distribution N (0, σ 2 ij ).Although the noises of each component are assumed to be Gaussian, the composition noise of each series by finite mixture are non-Gaussian and can display more flexibility in modeling.This gives us much freedom to fit some complicated data set.Model (2.1) assumes that the mixture probabilities are the same across the panel.It seems that this condition is a bit strict.In time series, to compensate the internal dependence of the data set, we often require a bit longer series than that in the independent case so that we can obtain good parameter estimates.However, in traditional panel time series data, each of the series contains only a few observations.For example, in the gray-sided voles data, each series contains only 31 observations.This sample size is also a bit small for statistical analysis even in the independent cases.To obtain accurate estimation, it is reasonable to assume that all these series have the same mixing probability so that we can use all the information we can obtain.When series length is a bit longer, we can remove this assumption.We will relax this restriction in Section 5.
Another advantage of the MARP model is that this model will widen the regions of the first order and second order stationarity.Wong and Li (2000) proved that the necessary and sufficient condition for the first order is that the roots z 1 , . . ., z p of the equation 1− 0 all lie inside the unit circle.The necessary and sufficient condition of second order stationarity is For details, see Wong and Li (2000).These results show that the mixture of non-stationary process and stationary process can be stationary. Let Here the subscript L in θ NL stands for "last" component so that we can avoid the confusion of the parameters in the component.Let π = (π 1 , . . ., π N ) T .Let Z jt = (Z 1jt , Z 2jt , . . ., Z Njt ) T , Z j = (Z T j1 , Z T j2 , . . ., Z T jT j ) T and Z = (Z T 1 , . . ., Z T M ) T .The vector Z jt contains the unobservable random variable Z ijt where Z ijt = 1 when at time t, in the jth series, the X jt comes from the ith component of the conditional distribution function and Z ijt = 0 otherwise.Then the parameters are naturally divided into two groups, θ and π with Z together.With the help of Z ijt , we can rewrite model (2.1) as (2.2).
where ε ijt is the white noise process corresponding to the ith component of the jth series.The {ε ijt } are also independent for i = 1, . . ., N, j = 1, . . ., M and all t.If we study equation (2.2) closely, we will find that the MARP model is actually a mixture of N Gaussian AR models for each of the series.The conditional mean of series j is Since the conditional mean depends on the past information of the time series, the shape of the conditional distributions will change from time to time and can be uni-modal or multi-modal.It can be proved that for a stationary mixture autoregressive time series X t , the marginal distribution of X t is also a mixture distribution.Due to this nature, the marginal distribution could be multi-modal if the conditional distribution is multi-modal.McLachlan and Peel (2000) showed some examples of uni-modal and multi-modal distributions in static situations and Wong and Li (2000) demonstrated the presence of multi-modality in dynamic situations.Note that for mixtures, the conditional mean could take the value at a local minimum point of the conditional mixture density and therefore it could be a misleading predictor of the future values.
Since Z is the vector of unobservable random variables, maximizing the loglikelihood function directly from the samples would be difficult.To overcome this difficulty, Dempster et al. (1977) suggested the EM algorithm which substitutes the unobservable variables by their expectations and then maximize the loglikelihood function.Using the expectation and maximization iteratively, one can obtain the estimates of the parameters when they converge.To use the EM algorithm, we should first find the log-likelihood function.The conditional loglikelihood of the jth series involving Z ijt is given by The global conditional log-likelihood function given Z for the entire panel is Since the first derivatives of the log-likelihood with respect to the interested parameters can be obtained straightforwardly, we skip them here.With the first derivatives, we can apply the EM algorithm.The iterative procedure of the EM algorithm can be divided into two steps -E-step and M-step respectively.

E-step
In this step, the parameter vector θ are assumed to be known.Then the missing data Z will be estimated by their conditional expectation respecting to θ and the observation vector X = (X T 1 , . . ., X T M ) T where X j = (X j1 , . . ., X jT j ) T .Let τ ijt denote the conditional expectation of Z ijt .Similar to Le et al. (1996) and Wong and Li (2000), the τ ijt can be expressed as where i = 1, . .., N , j = 1, . .., M , t = 1, . .., T j and ε ijt has the same meaning as before.

M-step
In this step, we replace Z ijt by its conditional expectation τ ijt .Therefore, π i can be estimated by the mean of .
We estimate the parameter vector θ by maximizing the global log-likelihood function l.This can be done by setting the first derivatives to zero and then solve the system of equations.The system can be simplified as follows.Let where ( φ011 , φ111 , . . ., φp,N−1,M ) is determined by another system of equations: ( where i = 1, . . ., N −1, j = 1, . . ., M, k = 1, . . ., p.For the N th component, since the φ ki are common to all the series, the estimation equations are a bit different from (2.5) -(2.7).By pooling all available information, we can write the equation system of φ0N , . .., φpN , σ2 N 1 , . .., σ2 NM as follows: Studying the system of equations carefully, we find that, by fixing j, to the ith component, equations (2.6) and (2.7) contain p+1 linear equations with p + 1 unknown variables φ 0ij , . .., φ pij .Hence, the parameters of the ith component in the jth series, i = 1, . . ., N − 1 can be solved separately from the other components.Let A ij be the matrix below: and let φij and b ij be Then equations (2.6) and (2.7) and the solution of the system can be written in matrix form as in (2.11) and (2.12) respectively.
and therefore σ 2 ij can be easily evaluated with equation (2.5).
The analysis of panel time series often encounter the difficulty that the number of observations in each time series is limited.Many researchers pooled the data from different series which shared a common structure, see for instance, Hjellvik and Tjφstheim (1999), Stenseth et al. (1999), Yao et al. (2000) and Fu et al. (2002a).They obtained more accurate results with the increase of the sample size.For the parameters of the N th component, they also have common parameters φ kN , k = 1, . . ., p. Hence, we pooled the observations X jt and τ Njt , j = 1, . . ., M to estimate φ kN , k = 0, . . ., p. Equations (2.8) -(2.10) form a (1 + p + M )equation system with 1 + p + M unknown variables.Since this equation system is nonlinear, a numerical method is needed to find the roots of the system.However, by substituting (2.8) to (2.9) and (2.10), the latter two equations contain only 1 + p unknown variables which will be easier to solve than directly solving the 1+p + M equations together since p is much smaller than M in most of the cases.
To obtain estimates of the parameters, we iterate the E-step and M-step until convergence.The performance of the EM algorithm will be shown in Section 3. Louis (1982) proposed the missing information principle by which the observed information matrix I can be computed by subtracting I m , the missing information matrix from I c , the complete information matrix. (2.13) Since the variance matrix of the estimates θ equals to the inverse of the observed information matrix I, equation (2.13) can be used to find the theoretical large sample variance of the estimates.All the entries in I c and I m are available from the authors.Model selection is an important problem in time series modeling.However, since the length of each series in the panel data is often very short, a high order model is not preferred.Setting the order p = 1 and the number of components N = 2 will be a wise choice.Hence, we will not pursue this problem.

Simulation Study
To investigate the performance of the proposed estimation procedure, we do a Monte Carlo simulation1 .As mentioned before, for a panel data set which the length of each series is short, a low order model is preferred.Hence, we simulated the data set from a twenty-series panel model.Each series is a two-component AR(1) mixture time series.The model is where Z 1jt = 1 with probability π 1 if X jt is from the first component of jth series.Otherwise, Z 1jt = 0.The noise terms ε 1jt and ε 2jt followed two zeromean normal distributions with variances σ 2 1j and σ 2 2j .The parameter values of φ 01j , φ 11j , σ 2 1j and σ 2 2j for j = 1, . . ., 20 are listed in Table 1 respectively.The mixture probability π 1 = 0.4 and π 2 = 1 − π 1 = 0.6.The parameters of the common structure φ 02 and φ 12 were 1.5 and -0.5 respectively.We did four experiments with different series length to study the finite sample properties of the proposed estimation method.The series lengths T j are around 1100, 100, 50, 30 in these four experiments respectively.Each experiment was replicated 1000 times.The results are listed in Table 2.
.5 -0.4 -0.2 φ 0,1,11 φ 0,1,12 φ 0,1,13 φ 0,1,14 φ 0,1,15 φ 0,1,16 φ 0,1,17 φ 0,1,18 φ 0,1,19 φ 0,1,20 0.2 0.4 0.5 0.6 0.8 1.0   From the results in Table 2, we can find that the estimation results are very good when series length was longer than 100.In the case that series length was shorter than 30, the estimate method was still applicable, though the standard errors were a bit large.In this extreme case, one can obtain acceptable results by applying our model.Therefore, the lower bound of the series length to apply our model and estimation method seems to be about 30.When better result is required, a series length of about 50 for each of the series is recommended.
Using the information matrix (2.13), we also calculated the theoretical standard error for the case that the length of each series is about 1100.The results for the φ s and π are listed in Table 3.The results for the σ2 s are similar and we choose not to list them.The results are in general close to the empirical standard errors in Table 2 suggesting the accuracy of (2.13) in large samples.Hence, our estimation procedure seems useful and can provide good results.

A Real Example
This data set consists of the number of yearly trapped gray-sided voles (Clethrionomys rufocanus) in ninety-one different stations in forested regions of northern Hokkaido, Japan.Hokkaido, the northen-most island of Japan displays similar bio-geographically feature as the neighboring Asian mainland more than the other Japanese islands (Tatewaki 1958;Kondo 1993).The average temperature in August and February is above 20 • C and below −7 • C. The amount of snow in the winter is largely correlated to the distance to the sea shore.The forests in which the voles were trapped were under the management of the Forestry Agency of the Japanese Government.Among them, 76 % is indigenous forests and 24 % is planted forest.The gray-sided vole is a forested habitat mammal which can be found in both types of forests.The observation period spanning from 1962 to 1992, over 31 years.This data set has been widely studied in mathematical ecology and panel data literature such as Bjornstard et al. (1996) and Stenseth et al. (1996).The ninety-one series can be divided into three groups with 16, 41 and 34 series in group one, two and three respectively.For details, see Hjellvik and Tjφstheim (1999) and Hjellvik et al. (2004).Figure 2 shows the series in group two after logarithm transformation.Hjellvik and Tjφstheim (1999) developed an autoregressive model with an intersection term to interpret the correlation among these panels.Hjellvik et al. (2004) developed a nonparametric estimation method to model the intercorrelation and a test to check the linearity assumption.Jin and Li (2005) assumed that the different series were correlated contemporaneously and proposed a semiparametric model to fit the data set.All these papers did not take the multi-modal phenomenon into account.Stenseth et al. (1996) tested the linearity of the data set and found that it was rejected at 5 % level.But they still employed linear models to fit the data.However, after the linear fitting, they found a clear indication of bimodality which in fact violate their underlying assumption.The histogram (Figure 1) of four typical series partially shows the bimodal characteristics.In many biological papers (for example, Bjornstard et al. (1996); Stenseth et al. (1996);Ota (1984) which studied the gray-sided voles, it has been pointed out that the population of the voles not only correlated with the previous population size but also correlated to the distance from the Sea of Japan.The latter factor affects the climate directly and the food of the gray-sided voles indirectly.The impact of these unobservable factors will randomly vary from year to year and they often followed different patterns.Hence, it seems appropriate to fit the data with MARP so that we can catch these features.Similar to Hjellvik and Tjφstheim (1999), we also used the series in group two which contained 41 series.The length of each series was 31.Among the 1271 observations, there were 174 (13.7%) zeroes.Hjellvik and Tjφstheim (1999) proposed the model (4.1)

Histogram of Voles1
to model the inter-correlation effect of this data set.Here φ k denotes the common autoregressive structure and η t is the possibly nonlinear effect influencing all the series at time t.In their analysis, they first treated η t as nuisance parameters and estimate them by the averages of X it over i and kept t fixed.They discarded the nonlinearity property since they believed that it was very weak.Clearly, the assumption that all the series follows the same autoregressive model is too restrictive and ignoring the nonlinearity may be inappropriate.Slightly loosening the restriction of the model would help us to find new implicit pattern inside the data.In fact, by ignoring η t , model (4.1) is a special case of model (2.1) or (2.2) if we set the probability π N = 1 − N −1 i=1 π i = 1.Although we have increased the number of parameters, we have gained much flexibility by relaxing the restrictive assumption of model (4.1).
Similar to Hjellvik and Tjφstheim (1999), for convenience, we added a value of one to each datum point and then took logarithm.After that, we fitted the data with model (2.2).In fact, we only used the first 30 observations in each series to fit our model.The last observation in all the series are used to check the coverage of our model.Table 4 shows the result.Since all these voles were found on the same island, it was reasonable to believe that they should follow some common pattern which were modeled by the common mixture component whereas the regional variations may be reflected by the individual component.
From the results in Table 4, we found that the autoregressive coefficients of the common structure are significantly different from those of the first component in each of these 41 series which suggested that equation (2.1) could be an appropriate model.Together with the result of π = 0.3015, we believed that the mixture phenomenon was quite significant which was consistent with the histogram plots (Figure 1) of these series.An interpretation of the result is that the data follow the common component about 70 % of the time whereas in the other 30 % time, the dependent structure follows a pattern which is dictated by the environment, possibly the food supplies.From our analysis, we found that among these 41 series, the individual components showed three different patterns.The twentyfirst, thirtieth and thirty-fifth series displayed negative autocorrelation in the individual components.These suggested that the regional factors of these three districts may have negative effect on the number of the voles.It is possible that the strict survival environments could not support too many voles.The other twenty-nine series showed positive autocorrelation in the individual components suggesting that the environment could be relatively good.The food might be sufficient and competition within species might be mild.The remaining nine series had insignificant autocorrelation with the previous observation.The φ 11j were around zero suggesting that the number of voles is a constant.To compare our model with the one of Hjellvik and Tjφstheim (1999), we checked the coverage of the one-step forecasting.Since our model was a bit complicated, we had to use simulation to find the empirical quantiles.For example, let X 1,30 be the thirtieth observation of the first series.From the results in Table 4, we knew that X 1,31 would follow a mixture distribution.It would come from the N (0.2316 + 0.2851X 1,30 , 0.7040) distribtyion with probability 0.3015 and come from the N (1.1513 + 0.6402X 1,30 , 0.3337) distribution with probability 0.6985.Note that all X jt have been log-transformed.Therefore, we can generate a random number U following a uniform (0,1) distribution.If U was smaller than 0.3015, we picked X1,31 from the first normal distribution as a simulated X 1,31 .Otherwise, we picked X1,31 from the second normal distribution.Repeating this 10,000 times, we obtained a sample following the mixture distribution.The empirical 95% predictive interval was obtained by (a, b) where a was the empirical 2.5% quantile and b was the empirical 97.5% quantile.With the same method, we obtained the empirical predictive intervals for all these 41 series.Comparing to the corresponding interval, we found that there were 3 last observations among these 41 series outside 95% one-step ahead predictive intervals, a coverage of about 7.3%.Considering the small sample size, this coverage was still very near the nominal 95%.There were 3 different model orders in Hjellvik and Tjφstheim (1999), one, two and three respectively.In general, higher order model would be more accurate than lower order model.Hence, comparing our model with their model whose order was three could be more cogent than comparing to the other two.By checking the coverage of their model, we found that there was no observation outside the 95% intervals.This indicated that their model leads to conservative predictive intervals.Another commonly used method for model comparison is BIC (Schwarz, 1978).This criterion adds heavier penalty than AIC (Akaike, 1974).Since the number of parameters are large, using BIC can be a challenge to our model.Although panel data are multivariate time series, since the residual of our model were hard to define due to the characteristics of the finite mixture, we used a simpler version of BIC which was the same as that in the univariate case, BIC = −2 log(maximum likelihood) + log(number of data points) ×(number of estimated parameters).
The values of the log-likelihood were -1082.10 and -1310.27for our model and that of Hjellvik and Tjφstheim (1999) respectively.The BICs were 2680.22 and 2715.40 respectively for these two different models.Although the difference between the two BICs are not very large, it is an evidence that our model better fits the data.

Extension
In this section, we relax the assumption that the mixture parameters π i are common to every series.For simplicity, we only show the two component and order one case.The more general cases can be extended straightforwardly.For notation consistency, we still use three indices as the subscripts of the parameters.The model is modified as where π 1j + π 2j = 1 for each j, j = 1, . . ., M. The conditional log-likelihood functions, for each j given Z j which has the same meaning as section 2, are Since the first derivatives of the log-likelihood function l can be obtained straightforwardly, we skip them here.We show the estimation method with the EM algorithm briefly since the method is almost the same as the one in section 2.

M-step
Suppose that the missing data are known.The estimates of the parameters θ can then be obtained by maximizing the log-likelihood function l.
The M-step equations become πij = T j t=p+1 τ ijt T j − p . (5.3) The estimation equations of the other parameters are the same as those in section 2.
Following the missing information principle Louis (1982), the information matrix equation (2.13) still holds.We skip the derivation of the entries in I c and I m since the method is just the same except for those entries involving π ij .To investigate the performance of model (5.1), we did a simulation study.The underlying model considered was (5.4) where Z 1jt = 1 with probability π 1j if X jt was from the first component of jth series.Otherwise, Z 1jt = 0.The parameter values of φ 0ij , φ 1ij , σ 1j , σ 2j and π 1j are listed in Table 5.The parameters of the common structure φ 02 and φ 12 were 1.5 and -0.5 respectively.It can be seen that the parameters are almost the same as model (3.1) except the parameters π 1j .We also did four experiments with different series length to study the finite sample properties of the proposed estimation method.The series length is also the same as before, say around 1100, 100, 50 and 30 respectively.Each experiment was replicated 1000 times.
The results are listed in Table 6.From the results, we find almost the same conclusion as Table 2.The difference is that the standard deviations of the φ j and σ 2 ij estimates become larger in Table 6 than those in Table 2.However, the bias of the estimates seems acceptable.We also applied model (5.1) to the gray-sided voles data.As before, we only used the first 30 observations in each series and left the last one to check the coverage of the forecasting value.The results are listed in Table 7.The results of series 8, 9 and 38 showed a very large π which were close to 1. Considering the short series length, we can accept some bad estimates.For series 8, the second component has a very large variance about 14.Note that the mixture probability of this series is 0.95, there were only two observations coming from the second component.We checked the two observations and found that they were exactly the two missing data points in this series.Hence, the variance of the second component was just an artifact from the Fortran program and would be meaningless.Other series showed a significant mixture phenomenon.By checking the coverage of the forecasting, we found that there were also 3 observations outside the 95% one-step ahead predictive interval, that is, a coverage about 7.3%.However, although the log-likelihood was -963.38 which were smaller than that of model (3.1) and the one in Hjellvik and Tφstheim (1999), the number of parameters was 207 which increased the BIC to 3399.52.Because of the large number of parameters, we have some reservations on fitting gray-sided voles data with model (5.1).

Conclusion
In this paper, we considered a mixture autoregressive panel (MARP) model which can capture the multi-modal phenomenon in some data sets.We suggested to estimate the MARP model by the EM algorithm.The simulation studies showed that we can obtain quite satisfying estimates if the series length of each series is 50 or longer.Even if the series length is around 30, the estimates still appeared to be acceptable.We illustrated the model by applying it to the graysided voles data.All these results show that the sample size requirement of our model is quite low and the performance of the model is satisfying so that it may be of some potential in applications.

Figure 1 :Figure 2 :
Figure 1: Histograms of four typical series in the gray-sided voles data.The solid lines stand for the kernel smooth of the data with a Gaussian kernel.

Table 1 :
True values of the parameters of model (3.1)

Table 2 :
Estimation results of the simulation of model (3.1) with different series length (1000 replications) 1000 < T j < 1100 80 < T j < 100 40 < T j < 50 20 < T j < 30 Only the results of series 1-4 and the common structure are shown in this table.Other estimates of the parameters can be obtained from the authors.est.means the average value of the estimates std.means standard deviation of the estimates.

Table 4 :
Estimation results of the gray-sided voles