1

CURRENT POPULATION SURVEY STATE-LEVEL 

VARIANCE ESTIMATION 

Richard Griffiths and Khandaker Mansur 

1

U.S. Bureau of the Census 

Abstract 

The problem of estimation for small area means and totals is well-documented in the literature. Not so, the problem 

of estimating variances for small areas. Yet, small area variance estimators suffer from the same general problem as 

small area estimators of means and totals – a small sample size on which to base the estimates. Current Population 

Survey state-level variance estimators are an example of variance estimation with relatively small sample sizes. In 

this paper we examine two modeling approaches to address this small area estimation problem. 

Key words: variance, Current Population Survey, small area estimation, modeling 

Introduction 

The Current Population Survey (CPS) sample design is a two-stage stratified, cluster design for 

each state and the District of Columbia. Within each state primary sampling units (PSUs), which 

are groups of counties, are stratified. A single PSU is selected into the sample in each stratum 

and a systematic sample of clusters of housing units is then drawn from each sampled PSU. 

Sampling is done independently in each state. 

There are two types of strata in the CPS sample design: self-representing (SR) and non-selfrepresenting (NSR). Each SR stratum contains a single PSU, which is selected into the sample 

with probability one. Each NSR stratum contains at least two PSUs, one of which is selected 

into the sample. The variances of CPS estimators thus have two components in NSR strata: a 

between-PSU component and a within-PSU component. In SR strata the estimators have only a 

within-PSU component of variance. 

The U.S. Census Bureau currently calculates monthly estimates of variances for CPS state labor 

force estimators. Both successive difference replication and modified half-sample replication 

methods are used to calculate these state-level variance estimates. (See Fay and Train, 1995, and 

U.S. Census Bureau, 2000.) The method of successive difference replication is used to estimate 

within-PSU variances in SR strata and the half-sample replication method is used to estimate 

total variances in NSR strata. These variance estimates suffer from two known problems: They 

are based on relatively small sample sizes and they are subject to a bias induced by the procedure 

of collapsing NSR strata to estimate between-PSU variances. In this paper we address the first 

of these problems by discussing methods for modeling the variance estimates to improve their 

precision. Griffiths and Mansur (forthcoming) relate a method for reducing the bias of the 

variance estimators. 

This paper reports the results of research and analysis undertaken by Census Bureau staff. It has undergone a 

Census Bureau review more limited in scope than that given to official Census Bureau publications. This report is 

released to inform interested parties of ongoing research and to encourage discussion of work in progress. 
The approach we take to modeling is colored by the purposes for which the state-level variance 

estimates are (or can be) used. One use is as vital signs of the CPS. Since the U.S. Census 

Bureau is responsible for the statistical methodology of the CPS, it is important for us to monitor 

the quality (the “health,” if you will) of the estimators that come from the CPS. The state-level 

variance estimates are important quantities in diagnosing the health of the CPS methodology – 

they are vital signs. A second use is for research purposes: accurate state-level variance 

estimates are important in designing changes to the CPS methodology. Another use of the 

variances is in models based on CPS data. In particular, state-level small area estimation models 

need good estimates of variance. Beyond that, small area models often require covariance 

estimates. (See, for example, Tiller, 1992.) 

From these uses we can see a need for modeling both variances and covariances. Our primary 

concern is in modeling the variances, as they are most useful to us as vital signs of the CPS. 

However, we recognize the need for estimating covariances – a need that will probably grow as 

small area estimation proliferates. The models we describe in this paper are, in one form or 

another, capable of modeling both variances and covariances. 

In the following sections, we describe a model introduced by Otto and Bell (1995) for modeling 

variance-covariance matrices for the Census Bureau’s Small Area Income and Poverty Estimates 

program state-level model. We then proceed to examine a generalized linear model (GLiM) fit 

by maximum partial likelihood estimation. This model explicitly models only variances, but by 

making some assumptions about the relationship between variances and covariances, could be 

used to estimate covariances. After introducing these models, we examine the results of fitting 

them to several years of monthly state-level variance and covariance estimates. 

The Otto/Bell Model 

Otto and Bell (1995) proposed a model for improving estimates of state-level variancecovariance matrices for March CPS income and poverty estimators. This model is based on the 

principle that the mean of the variance-covariance matrix estimator is a function of several 

components: state effects, characteristic estimates, and sample sizes. Furthermore, the model 

assumes a structural relationship between variance and covariance estimates based on the 

autocorrelated nature of CPS sampling errors. Otto and Bell (1995) assumed that the variancecovariance matrix estimator follows a Wishart distribution; inference about the mean of this 

distribution can be made either through maximum likelihood estimation or in a Bayesian 

framework by assuming a prior distribution on the model parameters. 

The Otto/Bell model assumes that ν sCs ~ Wishart(ν s , Vs ) , where Cs

 is the sample-based 

variance-covariance matrix estimator for state s, ν s is the degrees of freedom for Cs, and Vs is 

the mean variance-covariance matrix for state s. Cs is a MxM matrix, with M being the number 

of months for which we have estimated variances and covariances. We denote the (i,j)th element 

of Cs by Csij, where Csij is the covariance of the estimators from months i and j in state s. 

We assume Vs has the form 

V s (η ) = w s⋅ GVF s ⋅ R ⋅ GVF s⋅ w s 

~ 

(1)
where η is a vector of unknown parameters contained in the terms on the right-hand side of (1); 

~

ws is a diagonal matrix, with each diagonal element being the square root of the state effect for 

state s; GVFs is a diagonal matrix containing the square roots of the generalized variance 

function (GVF) estimates divided by the sample size for state s; and R is a matrix which 

accounts for the autocorrelated nature of the sampling errors. In this paper we assume the tth 

diagonal element of where a and b are GVF 

coefficients, xst is the estimated characteristic total in state s for the tth month, and nst is the state 

sample size in month t. This is the form of the official CPS GVF, divided by the state sample 

sizes. (See U.S. Census Bureau, 2000.) 

GVFs is a ⋅ x

2

st / nst + b⋅ xst / nst , t=1,2,...,M, 

The R matrix has the following form 

⎛ 1 Corr (e s1 ,e s2 ) Corr ( e s1 ,e s3 )  Corr (e ⎜

s1 ,e sM ) ˆ

÷ Corr (e s2 ,e s1 ) 1 Corr (e s2 ,e 

R s

) Corr ⎜ 3  (e s2 ,e sM )÷

= ⎜ ÷

⎜ ÷ ⎝ Corr (e sM ,e s1 ) Corr ( e sM ,e s2 ) Corr (e sM ,e s3 )  1  ̄

where est is the sampling error for the estimator from state s in month t. So, Corr(est,es,t-k ) is the 

lag k sampling error autocorrelation. We assume the sampling errors represent a stationary 

stochastic process. The process assumed then determines the form of the elements of R. As an 

example, if we assume the sampling errors follow an ARMA(1,1) process, then 

th where {εst }t=1,..., M is a white noise process, and the (i,j) element of R 

has the form: (See Vandaele, 1984, pp. 46-47.)

est = φ es , t −1 + εst − θ εs, t −1 , 

(1− φ θ )(φ − θ ) j−i 

2 φ −1

, i...j. 

1+ θ − 2φ θ 

We thus see that the Otto/Bell model assumes the mean of the variance estimator is the product 

of the state effect and the GVF divided by the sample size: E( C stt ) = w s (ax 

2

st / n st + bx st / nst ) ,

where ws is the state effect for state s. And it assumes that the mean of the covariance estimator 

Cstu, t...u, is ws ( ax 2/ n + bx / n )(ax 2

st st st st su / n su + bx su / n su ) ⋅ Corr (e st ,esu ).

Parameter Estimation 

The vector of unknown parameters in (1) may be written as 

η = (w 1 ,w 2 ,..., w 51,φ 1,φ 2 ,..., φ p ,θ 1,θ 2 ,..., θ q ,a ,b , df )

~ 

where the φ and θ parameters are from the ARMA process used to describe the sampling error 

autocorrelations, which will include seasonal terms, and df is a parameter used to estimate the 

degrees of freedom in each state ( ν s = df − hs , where hs is the number of strata in state s). In this 

paper, we examine using estimated monthly state-level variance-covariance matrices Cs to 

calculate maximum likelihood estimates of η and thus of the Vs. We note that Otto and Bell 

~

(1995) treated the ws as random effects in their work, because they had only five annual variance 

estimates from each state for model fitting. We treat the ws as fixed effects since we have more 

observations, owing to the fact that we fit the model using monthly estimates. 
Since we have assumed a Wishart distribution for ν s Cs , the likelihood function we work with is 

[det(ν C )] (ν s − M −1) ⎡ 1 −1  ̆ 

s s

exp − tr(V ν C ⎢  ̇ ⎣ 2 s s s )

 ̊

L(η |Cs ) = M 

~ 2 ν s M /2 π M ( M −1)/ 4 [det(V )] ν s /2

s ∏ Γ [(ν s − i + 1) / 2] 

i=1 

Since CPS sampling is done independently in each state, we assume the state-level variance 

estimators are independent and the likelihood function using data from all states is then given by 

where 

51 

L(η | C) = ∏ L (η | Cs) , 

~ s =1 ~ 

C=(C1,C2,...,C51). 

(2) 

Maximizing (2), or the log of (2), for η gives us η , the MLE. The MLE of Vs

 is then obtained 

~ ~ 

by substituting η into (1). 

~ 

Partial Likelihood Generalized Linear Model 

As an alternative to the Otto/Bell model for improving state-level variance estimates, we 

consider a generalized linear model (GLiM) fit through partial likelihood (PL) estimation. (See 

Kedem and Fokianos, forthcoming, and Fokianos and Kedem, 1998, for more on PL estimation 

in the context of GLiMs.) We examine this alternative model for two reasons: 

• The Otto/Bell model fits both variances and covariances. For some purposes (e.g., variances 

as CPS vital signs), we are more concerned with the estimation of variances than 

covariances. (In fact, in this paper, we emphasize variances over covariances.) For these 

purposes, a properly-specified model that fits only variances will be more efficient than a 

variance-covariance model. 

• This alternative model will help us evaluate the fit of the Otto/Bell model for variances. 

The basis of this method is a decomposition of the joint density of the Cstt and xst. We may write 

the joint probability density function for state s as 

f s(C s11 ,C s22 ,..., C sMM , x s1 , x s2 ,..., x sM ,A s ) = 

M M 

f s (xs1|As )∏ f s (xst |Cs11, xs1,..., Cs , t − 1,t −1, xs , t − 1, As )∏ fs (Cstt |Cs11,..., Cs , t −1,t −1, xs1,... xst , As )

t= 2 t =1 

(3)

where s is the fixed auxiliary information for state s. Here s=(ns1,...,nsM). Rather than basing 

inference on the likelihood, we will base it on the partial likelihood. The last product in (3) is the 

partial likelihood. 

A A

Since CPS sampling is done independently in each state, we take the partial likelihood over all 

states to be 

∏

51 

∏ 

M 

f

s (Cstt |Cs11,..., Cs , t −1,t −1

, xs1,..., xst ,As) 

s=1 t=1 

(4) 
To develop the model for state-level variances, we assume that fs(Cstt|Cs11,...,Cs,t-1,t-1,xs1,...xst,As) 

is a gamma density, s=1, ..., 51; that Vstt is the conditional mean of the state-level variance 

estimator for month t: V stt = E( C stt |C s11 ,..., C s, t −1,t−1 , x s1 ,..., xst ) ; and that

g( V 2

stt ) = w s + β 1 C s, t −1,t−1 + β 2 x st / n st + β 3 x st / nst , 

for some link function g. We include the most recent lagged value Cs,t-1,t-1 in the formulation 

based on evidence that state-level variance estimates follow an AR(1) process. (See Mansur and 

Griffiths, 2001.) Since Vstt is a parameter in the density, (4) is a function of 

β = (w 1 ,.., w 51 ,β 1,β 2, β 3 ) 

~

. Specifying a form for g allows us to estimate β , and thus Vstt, by

~ 

maximizing (4) for β . This estimator is the maximum partial likelihood estimator (MPLE) of 

~ 

β . (See Wong, 1986, and Fokianos and Kedem, 1998, for discussions of the asymptotic 

~ 

properties of the MPLE.) 

Choosing the link function as g(Vstt) = ln(Vstt) gives a model with multiplicative state effects: 

V = e

w

2

s e 

β1Cs , t − 1,t−1 +β2 xst /nst +β3 xst /nst

stt . 

Choosing the link function as g(Vstt) = Vstt (the identity link) gives a model with the mean as a 

linear function of the covariates: 

V stt = w s + β 1 C s, t − + β 2

1,t−1 2 x st / n st + β 3 x st / n st .

Both forms of link function exhibit commonalities with the specification of the mean variance in 

the Otto/Bell model and thus will be used in comparing the PL GLiM to the Otto/Bell model in 

this paper. Note that the PL GLiM gives us a model very similar to a GVF and provides a 

theoretical justification for regressing on random covariates and past variance estimates. 

Results from Fitting the Models 

In this section we discuss results obtained by fitting the models to state-level total variance and 

covariance estimates. We used estimated variances and covariances for the (uncomposited) 

estimates of number of people unemployed from January 1996 to September 2000 for this fitting. 

We made no attempt to account for the bias due to collapsing NSR strata. 

We first examined the fit of several versions of both the Otto/Bell model and PL GLiMs. After 

determining which of these models were most appropriate to use in a further evaluation, we 

compared their fits for the estimated variances. 

Fit of the Otto/Bell Model 

We examined the fit of the full model (1) and several reduced models. The R matrix we used in 

fitting these models is for an ARMA(1,1)x(0,1)12 process. This is the form that corresponds to 

sampling error autocorrelation patterns for estimated number of people unemployed. (See 

Griffiths and Mansur, 2000.) We looked at the following reduced models: 

• Model (1) with Θ = 0 (i.e. no seasonal MA term) 12 

• Model (1) with ws = 1 for all s (i.e. no state effects) 
 

Table 1 Results of Fitting Full and Reduced Otto/Bell Models 

Model Number of Parameters 

Estimated2 

AIC 

Full Model 48 241,322 

Model with no seasonal MA term 47 241,976 

Model with no state effects 6 305,254 

Table 1 gives the AIC values for the fit of each of these models, along with the number of 

parameters estimated for each model. The AIC is calculated as , where k is

the number of parameters in η , which depends on the version of the model being fit, and S 

~ 

indexes the states used in the fitting (see footnote). 

−2∑ ln L(η | Cs ) + 2k 

s∈S ~ 

From this table we conclude that the state effects and seasonal MA term are important in the 

model. We thus used the full version of model (1) in the subsequent evaluation. 

Fit of the Partial Likelihood GLiMs 

We performed a similar analysis on the PL GLiMs with log link and identity link functions. 

Tables 2 and 3 give the AIC values and number of parameters estimated for each model. For 

these models we calculated the AIC as − 2 ln( PL) + 2k , where PL is the partial likelihood given 

by (4) conditioned on an observation at time t=0 and k is the number of parameters estimated for 

the model being fit; thus, 

M 

AIC = − 2∑ ∑ ln f s (C stt |C s00 ,..., C s, t −1 ,t − 1, xs 1 ,..., x st , As ) + 2k

s∈S t =1 

where Cs00 is the estimated variance from the month prior to month t=1. From the results shown 

in these tables, we conclude that the full models are the appropriate models to use in evaluating 

the fits of the PL GLiMs. 

Table 2 Fit of Full and Reduced Versions of PL GLiM with Identity Link 

Model with identity link Number of parameters 

Estimated1 

AIC 

Full Model 45 25,709 

Model with no lagged value of Cstt 44 25,736 

Model with no state effects 3 26,281 

Table 3 Fit of Full and Reduced Versions of PL GLiM with Log Link 

Model with log link Number of parameters 

Estimated1 

AIC 

Full Model 45 25,951 

Model with no lagged value of Cstt 

Model with no state effects

44 

3 

25,959 

33,755 

2

 The number of parameters excludes state effects for states with no NSR PSUs and Hawaii. These states were not 

included in the model fitting. 
Comparison of the Models for Variances 

Comparing the fit of the models for improving variance estimates is in some sense unfair since 

the Otto/Bell model was designed to model both variances and covariances and the PL GLiMs 

were designed to model only variances. On the other hand, the form of the GLiMs studied in this 

paper was restricted to make them somewhat comparable to the Otto/Bell model, at least in terms 

of link functions used. However, the comparison is made to help determine if we need different 

models for the different purposes mentioned in the opening section of this paper, or if one model 

will work for all our purposes. Below, we first compare the models on deviance. We then 

examine the models in terms of the similarity of their fits. To do this we look at the degree of 

smoothing each model provides and at graphs of the observed and modeled variances. 

Deviance 

The deviance is a measure of the discrepancy between the modeled and observed variances – a 

goodness-of-fit statistic. To determine the appropriate form of the deviance to use in comparing 

the models, we note that we assumed a gamma distribution for the Cstt under the PL GLiMs. 

Under the Otto/Bell model where 2 ν s 

s C s ~ Wishart ( ν s , V s ) ; thus, ν sC stt / V stt ~ χ (ν s ) = gamma( , 2) ,

2 

ν 

V is the tth

stt diagonal element of Vs

. So, we examined the deviance of the model fits under the 

assumption of a gamma distribution on the Cstt. 

The model deviance for a gamma distribution may be written as 

M 

2∑ [− ln(C /V  

stt stt ) + (Cstt − Vstt )/ V stt ], 

t =1 

where V 

stt is the modeled variance for state s, month t. (See McCullagh and Nelder, 1983.) We 

calculated the deviance of each of the models for all states. We found that the overall deviance 

for the Otto/Bell model (230.3) was similar to that of the PL GLiM with identity link (233.8), 

while that of the PL GLiM with log link (259.1) was quite a bit larger. A detailed look at the 

deviances by state also indicated that the Otto/Bell model and PL GLiM with identity link had 

generally similar deviances, while the PL GLiM with log link tended to have consistently larger 

deviances. Based on the deviance criteria, then, we might conclude that the Otto/Bell model and 

the PL GLiM with identity link fit the observed state-level variance data about equally well. 

Similarity of the Fitted Variances 

Graphs of the time series of observed state-level variances and modeled variances under each of 

the models are given in Figure 1 for several states. Two things stand out in these graphs: the 

modeled variances are similar under all three models and they generally describe smoother time 

series than do the observed variances. The similarity of modeled variances, especially those of 

the Otto/Bell model and the PL GLiM with identity link, is a result of the commonalities among 

the models, notably in the form of the hypothesized means of the variance estimators. 

As for the smoothing, under the assumption that the variance estimates are somewhat unstable 

(i.e., the small sample size problem), we would like to see time series of modeled variances that 
Table 4 Degree of Smoothing 

Characteristic of the 

Distribution 

Otto/Bell Model Partial Likelihood GLiM 

with identity link 

Partial Likelihood 

GLiM with log link 

Mean .466 .477 .389 

Median .434 .473 .334 

Standard Deviation .096 .136 .215 

are less oscillatory than those of observed variances. To measure the smoothness of each 

variance time series, we looked at the total variation in the observed and modeled variance time 

series. We defined the total variation in a time series 

M 

{Vt}t=1,...,M as ∑ V t − V t −1 . 

t =2 

We used the ratio of the total variation in the modeled variances to the total variation in the 

observed variances to assess the degree of smoothing. We calculated this ratio for each state and 

each model. Table 4 contains some characteristics of the distribution of this ratio over all states 

for the models. 

Before proceeding, we note the informal nature of the total variation in assessing the smoothness 

of the fits. We have no absolute criterion for the degree of smoothing that is best. It is certainly 

possible to smooth the estimates too much. In other words, there could well be a good deal of 

oscillation in a time series of true variances (e.g., seasonality) and smoothing the variances too 

much might mask this true oscillation. However, we feel this measure of smoothness provides a 

nice descriptive tool for understanding the graphs in Figure 1, plus it gives us another way of 

showing the similarities among the modeled variances. 

Table 4 shows that the degree of smoothing attained by the PL GLiM with log link was, on 

average, greater than that of the other two models, as well as more variable over the states. 

Overall, the smoothing for the PL GLiM with identity link was similar to that of the Otto/Bell 

model, though a little more variable over the states. This is another indication of the similarity of 

the fitted variances from the Otto/Bell and PL GLiM with identity link models. 

Discussion 

In this paper we have examined models for improving CPS state-level variance and covariance 

estimates, though we have been primarily concerned with the variance estimates. We have seen 

that fitting the Otto/Bell model to variance-covariance matrices and a PL GLiM to variance 

estimates resulted in substantially similar modeled variances. 

While the Otto/Bell model and the PL GLiM are aimed at somewhat different goals, their 

structures, along with some of the assumptions we’ve made, will allow us to substitute one for 

the other. We have seen that the Otto/Bell model when examined solely for its fit to variance 

estimates performed at least as well as the PL GLiMs examined in this paper. It would thus 

seem that, unless we can specify a much better form for the PL GLiM, the Otto/Bell model might 

be preferred for modeling both variances and covariances, purely on the basis of quality of fit. 

With this in mind, though, we also note that given a set of modeled variances from the Otto/Bell 

model, the modeled covariances will be completely determined by the estimated ARMA 
parameters in the R matrix. Thus, if we assume the same structural relationship between 

variance and covariance estimates as the Otto/Bell model, we can use the PL GLiM variances to 

calculate covariance estimates. Since we have seen that the Otto/Bell model and PL GLiM 

produced similar variance estimates, applying the same estimated ARMA parameters to the PL 

GLiM variance estimates would give modeled covariances similar to those of the Otto/Bell 

model. Note, though, that without first fitting the Otto/Bell model, the ARMA parameters would 

have to be estimated outside of the PL GLiM model fitting. They would have to be based on an 

analysis of sampling error autocorrelations similar to that of Griffiths and Mansur (2000). 

Finally, we note that we need to do more work to determine an appropriate link function for the 

PL GLiM. The links studied in this paper were used more for comparability with the Otto/Bell 

model than for quality of fit. However, due to the relative computational simplicity of 

implementing the PL GLiM over the Otto/Bell model, we believe that a form of PL GLiM 

should be used to produce modeled CPS state-level variance and covariance estimates. 

Acknowledgments 

The authors would like to thank Bill Bell and Harland Shoemaker of the Census Bureau for their 

helpful comments. 

References 

Fay, R.E. and G.F. Train (1995), “Aspects of Survey and Model-based Postcensal Estimation of 

Income and Poverty Characteristics for States and Counties,” Proceedings of the Section 

on Government Statistics, American Statistical Association, 154-159. 

Fokianos, K. and B. Kedem (1998), “Prediction and Classification of Non-stationary Categorical 

Time Series,” Journal of Multivariate Analysis, 277-296. 

Griffiths, R. and K. Mansur (2000), “Preliminary Analysis of State Variance Data: Sampling 

Error Autocorrelations (VAR90-36),” Internal U.S. Census Bureau memorandum. 

Griffiths, R. and K. Mansur (forthcoming), “The CPS State Variance Estimation Story (VAR90-

38),” Internal U.S. Census Bureau memorandum. 

Kedem, B. and K. Fokianos (forthcoming), Regression Models for Time Series Analysis, John 

Wiley & Sons, Inc. 

Mansur, K. and R. Griffiths (2001), “Analysis of the Current Population Survey State Variance 

Estimates,” paper presented at the 2001 Joint Statistical Meetings of the American 

Statistical Association, Section on Survey Research Methods. 

McCullagh, P. and J.A. Nelder (1983), Generalized Linear Models, Chapman and Hall. 

Otto, M.C. and W.R. Bell (1995), “Sampling Error Modelling of Poverty and Income Statistics

for States,” Proceedings of the Section on Government Statistics, American Statistical 

Association, 160-165. 

Tiller, R.B. (1992), “Time Series Modeling of Sample Survey Data from the U.S. Current 

Population Survey,” Journal of Official Statistics, 149-166. 

U.S. Census Bureau, Bureau of Labor Statistics (2000), Current Population Survey: Design and 

Methodology, Technical Paper 63,Washington, DC. 

Vandaele, W. (1983), Applied Time Series and Box-Jenkins Models, Academic Press, Inc. 

Wong, W.H. (1986), “Theory of Partial Likelihood,” The Annals of Statistics, 88-123. 
Observed 

4.50E+08 

Otto/Bell PL GLIM ID 

Iowa 

PL GLIM LOG 

9.00E+08 

South Carolina 

3.00E+08 6.00E+08 

1.50E+08 3.00E+08 

0.00E+00 0.00E+00 

Figure 1 Time Series of Observed and Modeled Variances 

5.00E+09 

1.20E+08 Texas Utah 

4.00E+09 

3.00E+09 7.00E+07 

2.00E+09 

2.00E+07 1.00E+09 Jan-96 Jan-96 Jan-97 Jan-97 Jan-98 Jan-98 Jan-99 Jan-99 Jan-00 Jan-00 

Jan-96 Jan-96Jan-97 Jan-97Jan-98Jan-99

Jan-99

Jan-00 

Jan-98Jan-00