HITOTSUBASHI

UNIVERSITY

Center for Financial Engineering Education (CFEE)CFEE Discussion Paper Series No. 2008-2

Bayesian Estimation of Unknown Regression Error Heteroscedasticity

Hiroaki Chigira and Tsunemasa Shiba

Center for Financial Engineering Education (CFEE)
Department of Economics, Hitotsubashi University
2-1 Naka, Kunitachi Tokyo, 186-8601 Japan
http://www.econ.hit-u.ac.jp/~finmodel

Bayesian Estimation of Unknown
Regression Error Heteroscedasticity
Hiroaki Chigira∗and Tsunemasa Shiba†
March 11, 2009

Abstract
We propose a Bayesian procedure to estimate heteroscedastic variances of the regression error term ω, when the form of heteroscedasticity is unknown. The prior information on ω is elicited from the wellknown Eicker–White Heteroscedasticity Consistent Variance–CovarianceMatrix Estimator. Markov Chain Monte Carlo algorithm is used tosimulate posterior pdf’s of the unknown elements of ω. In addition to
the numerical examples, we present an empirical investigation of thestock prices of Japanese pharmaceutical and biomedical companies todemonstrate usefulness of the proposed method.

key words
Eicker–White HCCM, orthogonal regressors, informative prior pdf’s,
MCMC, stock return variance
JEL Classiﬁcation
C11, C13∗

Department of Economics, Tohoku University, 27-1, Kawauchi, Aoba-ku, Sendai-shi
980-8576 Japan. e-mail: hchigira@econ.tohoku.ac.jp†
Department of Economics, Hitotsubashi University, Kunitachi-shi, Tokyo 186-8601
Japan. e-mail: tshiba@stat.hit-u.ac.jp Tel & fax: +81-42-580-8786

1

1

Introduction

We propose a Bayesian procedure to estimate heteroscedastic variances of
the regression error term, when the form of heteroscedasticity is unknown.
As pointed out by Amemiya (1985, p.199), the crucial ω vector1 cannot be
consistently estimated because as the number of parameters increases, thesample size also increases at the same rate, leading to the lack of identiﬁability of ω. In asymptotics framework, Eicker (1963) and White (1980)
independently developed a well-known consistent variance-covariance matrix estimator (“HCCM” hereafter) for the OLS regression coeﬃcient estimator. The methodology we propose in this paper is a Bayesian that uses
information obtained from the HCCM, in terms of a candidate density of a
Metropolis-Hastings (“M-H” hereafter) algorithm in Markov Chain Monte
Carlo simulation. The lack of identiﬁability of ω poses no problem. For one
thing, as in Amemiya (1985) we use an orthogonal regression that circumvents possible underidentiﬁability of ω, and we shall explain this method in
detail later2 . Second, we impose a prior on ω so that the vector becomes
identiﬁable in a Bayesian context.
The trend in the HCCM literature seems to be how to improve the ﬁnite
sample performance of tests of the linear restriction(s) on the coeﬃcient
vector, e.g., Long and Ervin(2000) and Godfrey (2006), among others. We
note that our focus in this paper is in the direct estimation of the elements
ω vector. There are papers that deal with statistical inferences of regression
coeﬃcients, when the scedastic function of the error term is unconstrained.
Robinson (1987), for example, assumes it to be a function of regressors, and1
The ! vector has in its elements, all the normalized diagonal elements of variancecovariance matrix of the regression error term. The normalization rule for the matrix is
given just below equation (1).
2
Amemiya op cit proposes and uses an orthogonal regression to obtain a better performing GLS.

2

derives an GLS estimator that is more eﬃcient than the existing ones. Our
Bayesian estimation of heteroscedasticity should sharpen posterior densityof regression coeﬃcient vector β and/or lead to a better predictive density.
It may also lead to more eﬃcient estimator of β in terms of asymptotic
theory framework as well.
We need to discuss the direct estimation of the ω vector. In ﬁnancial
returns, ω is nothing but the volatility. In order to access an option pricing,what we need to do ﬁrst is to come up with a reasonable estimate of volatility.
Our estimation of ω needs no parametric model for the volatility process
such as the GARCH model, since we use information obtained from the
HCCM estimation, in our MCMC simulation. If we wish to estimate avolatility process in time series data nonparametrically, what we usually do
is to calculate a historical volatility series. But this is just a descriptive
statistic without a theoretical background. Moreover, when it comes tocross section data, historical volatility calculation breaks down for obvious
reasons. Our Bayesian method, on the other hand, should provide a gooddeal of theoretical support for cross sectional data.
Our strategy to estimate the ω vector is Bayesian. After assuming ausual prior density for the parameters in the regression model, we are able towrite down a joint posterior density. The usual parameters such as regression
coeﬃcients may be easily simulated using the Gibbs sampler scheme. It is in
the simulation of the elements of ω that we use the HCCM. We use resultsfrom HCCM to form the candidate density in the M-H algorithm.
The rest of this paper is organized as follows. In section 2, we set our
regression model. Prior pdf’s are assumed here, and the joint posterior pdf
is derived. Section 3 starts out with our Bayesian MCMC calculation by a3

Gibbs sampler. We propose to use the Eicker–White result to simulate ω bya M-H scheme. Our numerical illustration and an empirical investigation of
the stock prices of Japanese pharmaceutical and biomedical companies, aregiven in section 4. Section 5 concludes.

2

The Model and the Joint Posterior pdf

Suppose a heteroscedastic regression model,
y = Xβ + ε,(1)where X ∼ n × K is a purely exogenous variable matrix, ε ∼ N(0, σ 2 Ω),Ω = diag(ω) = diag(ω1 , . . . , ωn ), ω ⊂ n+ , tr(Ω) =ni=1 ωi = n, σ

2 > 0,and n is the sample size3 . Let the column vector θ be all the parameters in
the model, i.e., θ  = (β  σ 2 ω  ) ∼ 1 × (K + 1 + n). The likelihood function,
then, becomes−1/2 −np(y|θ, X) ∝ |Ω|

σ





1
exp − 2 (y − Xβ) Ω−1 (y − Xβ) .
2σ(2)Assuming that the components of θ are independent a priori, and specifically the prior of β is diﬀuse, we obtain the following prior pdf under thetr(Ω) = n assumption:p(θ) ∝ |σ 2 Ω|−1/2 1[ω,ω ] (ω),(3)

where “1[ω,ω ] (ω)” designates an indicator function such that it takes on thevalue one if an element of the ω vector is contained in the range [ω, ω], andzero otherwise.
We shall discuss how we obtained the upper bound ω and the lower
bound ω. Our scaling assumption of ω is given by3ni=1 ωi = n. This means“tr(Ω) = n” is a technical assumption that is often emplyed for identiﬁability purpose
in regression models with heteroscedastic error term. See, e.g., Greene (2008, p.158).

4

that the sample mean of ω is one, i.e., ω̄ = 1. Consider a rather extreme
case that ten percent of ωi ’s have ten times the magnitude of ω̄ = 1, i.e.,
10% of elements of ω is ωi = 10. This leaves the rest of ninety percentof ωj ’s, zeros only, i.e., ωj = 0 for j = i, where the i belong to the ﬁrst
ten percent of the entire n. Hence, along with the assumption ωi > 0 for
i = 1, . . . , n, the fact that ω̄ = 1 encourages us to set a certain range foreach ωi that depends on the sample size, n. In the rest of this paper, oursample size in simulation and also empirical investigation, is ﬁfty. We may,therefore, safely set (ω, ω) = (0.025, 12.5), say. What does this imply?
Consider an example of the ﬁrst four observations taking up the bulk of theentire volatility. That is4

ωi = 46 while

i=1

50


ωi = 4. If this is the case, ouri=5ω and ω may be given by4= 0.087 > 0.025 = ω
46

and

46
= 11.5 < 12.5 = ω.
4

The above implies that even in this extreme case, our lower bound, ω, is far
smaller than the average of the latter forty six ωi ’s, while our upper bound,ω, far exceeds the average of the ﬁrst four ωi ’s. We emphasize that this data
dependent prior on the bounds of ωi for i = 1, . . . , n, are relatively easy toset, since the only information needed to form it is n. We also emphasizethat our experiments show insensitivity of MCMC draws to the bounds,
in the sense that the results are robust to the values of approximately ten
percent changes in the values of (ω, ω). Our prior pdf for θ in this paper isgiven byp(θ) ∝ σ −1 |Ω|−1/2 1[0.025, 12.5] (ω),(4)

and this replaces equation (3). Note that this speciﬁcation of (ω, ω) amounts
to(ω, ω) = (0.025/50, 12.5/50) = (0.0005, 0.25)
5

that is ω = 0.05 percent of n, while ω = 25 percent of n. This is the ruleof thumb that we recommend. When implementing this “rule of thumb,”around ten percent of each bound, should be explored.
Our joint posterior pdf is proportional to the product of equation (2)times (4):p(θ|y, X) ∝ |Ω|−1 σ −n−1 exp −1(y − Xβ) Ω−1 (y − Xβ) 1[0.025, 12.5] (ω)
2σ 2
(5)

3

Posterior Simulation: MCMC

We shall simulate β and σ using a Gibbs sampler, and then combine a
Metropolis-Hastings (“MH” hereafter) algorithm to simulate the entire θ 4 .

3.1

Simulating βLet θ −˛ represent a vector θ without β. Then the Gibbs sampler for β
starts out with
(0)

(0)

β (1) |θ −˛ ∼ p(β|θ −˛ , y, X),
where the superscript r in β (r) denotes the rth simulated β (thus r = 0 is
the initinal value) for r = 1, . . . , R, and R is the number of replications in
the MCMC. The conditional probability in the right hand side represents
a fully conditional posterior pdf of β, where a notation “A ∼ p(·)” is used
to indicate that a value A is generated from p(·) pdf. We use the OLS
estimated residual vector e = y − X β̂ to obtain s2 = e e/(n − K), and thisis used as the initial value of σ 2 . Next, note that it is pretty reasonable toset the initial values of ω, ω (0) , as ω (0) = ιn , where ιn ∼ n × 1 is a vectorof all one’s. This comes from the scaling assumption of ω, i.e., ω̄ = 1. As a4

Following several sections, in particular the Gibbs sampler subsection, may seem tobe a textbook content. But the notations introduced in these subsections are needed in
the subsection where we outline simulation of ! (Section 3.3).

6

(0)

result, our initial value vector becomes θ −˛ = (s2 ιn ) . Notice that the fullconditional distribution of β given θ −˛ isβ|(θ −˛ , y, X) ∼ N β̂, σ −2 (X  Ω−1 X)−1 .
(0)Then, using θ −˛ , we may generate β (1) from a multivariate normal conditional density:

β (1) ∼ N β̂, s2 (X  X)−1 ,where it turns out that the diagonal elements of the initial value of Ω issimply diag(ιn ) = I n . In r ≥ 2 iterations, above β (1) generation scheme isreplaced by


2



−1

β (r) ∼ N β (r−1) , σ (r−1) (X  Ω(r−1) X)−1 ,(6)

however.

3.2

Simulating σ

We, next, turn to generation of σ. The conditional posterior pdf of σ becomesp(σ|θ −ﬀ2 , y, X) ∝ σ

−n−1

exp −

ε Ω−1 ε 2
/σ ,
2

(7)


where ε = y−Xβ. This is the kernel of an inverse gamma pdf, IG(0)n ε Ω−1 ε
,
2
2



The initial value of θ −ﬀ2 becomes θ −ﬀ2 = (β (1) ιn ) , hence σ (1) is generated
as follows:



σ

5

(1)

(0)
|θ −ﬀ2 ∼ IG−(α+1)

n ε(1) ε(1)
,
2
2

,

We denote f (x|α, β) ∝ xexp(−β/x)1[0,+∞) (x) distribution as IG(α, β), inverse gamma with parameters α and β. Such x has E(x) = β/(α − 1) for α > 1,
and Var(x) = β 2 /((α − 1)2 (α − 2)) for α > 2. See, e.g., Robert and Casella (2004,
p.582). Our fully conditional σ pdf is given in equation (7), however, has σ 2 in the
denominator in the exponential function. To match the inverse gamma expression to
the above, make a transformation of random variables from x to σ and we obtain−(2α+1)exp(−β/σ 2 )1[0,+∞) (σ). IG(·) expression in the text, thus, follows.
f (σ|α, β) ∝ σ

75.

where ε(1) = y − Xβ (1) . For r ≥ 2 simulations, the above is replaced by(r−1)
σ (r) |θ −ﬀ2 ∼ IG−1n ε(r) Ω(r−1) ε(r)
,
2
2

.

(8)

Simulating ω Using HCCM Information

3.3
3.3.1

Conditional pdf of ωi

Let us now derive conditional distribution, ωi |θ −ωi , y, X. It is easy to seethat the conditional distribution of ω is in the form,−1p(ω|θ −! , y, X) ∝ |Ω|

1
exp − 2 ε Ω−1 ε 1[0.025, 12.5] (ω).
2σInside the exponential function in the above, ε Ω−1 ε becomes ε Ω−1 ε =
n

t=1



εi
√
ωi

2


√
= ε∗ ε∗ , where ε∗ = (ε∗1 , . . . , ε∗n ) ∼ n × 1 and ε∗i = εi / ωi ,nsince Ω = diag(ω). Note also |Ω|−1 =

t=1

ωi ,



p(ωi |θ −ωi , y, X) ∝ ωi−1 exp−

ωi−1 , we then have for a particularε2i/ωi 1[0.025, 12.5] (ωi ),2σ 2where εi = yi − xi β and xi is the ith row of X. The above is not a tractablepdf for two reasons, however. First, if we did not have a bounds constraint,∞[ω, ω], this function cannot be integrated, i.e.,0ωi−1 exp−

ε2i/ωi dωi2σ 2is not ﬁnite. Secondly, even if we let the ﬁrst parameter of the inverse
gamma pdf to be some positive value no > 16 , it would be awfully diﬃcult
to simulate a truncated inverse gamma random variable that has parameters,
IG(no ,

ε2i), in a Gibbs scheme, i.e., always accepting whatever is generated.
2σ 2

3.3.2

HCCM Information to be used in M-H

Let us, for a moment, depart from the problem of generating ωi , and see
what the sampling theory has to tell us about the heteroscedastic ε. An
6
We note that setting such no > 0, when it is not, is quite arbitrary. To deal with thispossibly unidentiﬁed conditional posterior pdf, we propose to set bounds on the ! vector,
and employ an M-H scheme, as discussed in section 3.3.3, below. Instead of this avenueof approach, we might assume a prior on the ﬁrst parameter in the inserve gamma pdf,
to turn our model to a hierarchical Bayes model.

8

OLS estimated residual vector, e = (e1 , . . . , en ) , is already deﬁned forequation (1).

Eicker (1963) and White (1980)’s result may be used toform a sandwitch variance-covariance matrix estimator of β̂, HCCM, that
is consistent:

1 
XX
n

−1



Ĥ

1 
XX
n

−1

n
1
, where Ĥ =xi xi e2i andn i=1

X  = (x1 , . . . , xn ) ∼ K × n. Under the usual assumptions such as a full1
rank matrix lim X  X ≡ Q < ∞, it is well known thatnpĤ −→ H = lim

n
σ2 
σ2 
X ΩX = lim
xi xi ωi
n
n i=1

(9)

holds. We may, hence, postulate a multivariate regression model of thefollowing:Ĥ = H + V ,(10)where V ∼ K × K is an error term matrix that is constrained to make both
H and Ĥ > 0.

7

After some manipulations (see the appendix to thispaper), the above becomesĥ = Xn σ 2 ω + v,

(11)

where ĥ = vech(Ĥ) ∼ K  × 1, Xn = [vech(x1 x1 ) . . . vech(xn xn )] ∼ K  × n,
ω = (ω1 . . . ωn ) ∼ n × 1, and ω ⊂ n+ , v = vech(V ) ∼ K  × 1, and K  ≡12 K(K + 1). This is a regression of K

dimension ĥ variable on K  × n

explanatory variable matrix Xn . Hence the degrees of freedom, here, is df =K  −n, and under the usual situations such “df” could even be negative, i.e.,the ĥ regression is unestimable. We have now found that the ω parametersin equation (11) are identiﬁable only when the β vector in equation (1) is
identiﬁable. To restate,
the ω vector is identiﬁable if the inequalities K  > n and n > K
simultaneously hold.
7

We use the notation “A > 0” to denote that the matrix A is positive deﬁnite.

9In a typical case like n = 50 and K = 3, ω is unidentiﬁable since K  = 6,
however, take a case like n = 20 and K = 10 then K  = 55 which makes ω tobe identiﬁable. But this is, indeed, a rare case. ω may become identiﬁable
if the column dimension of X increases, then, the “sample size” in equation(11), i.e., K  , should also increase. What we need to do is to increase K 
without altering e, estimated regression residual vector from equation (1),that constitutes a part of dependent variable in equation (11).
Amemiya (1983) proposed “partially generalized least squares” estimation method to improve upon eﬃciency over OLSE. We may use his idea tocircumvent the problem that we stated in the previous paragraph. Let the
regressor matrix in equation (1) to be
Z = (X W ) ∼ n × (K + KW ),(12)

where we require that W is orthogonal to all the variables in equation(1), i.e., W  (yX) = 0 ∼ KW × (K + 1) 8 . It is easy to see that b

and e from equation (1) remain the same even after we regress y on Z
instead of X alone. We now note that the the identiﬁability condition
given in the previous paragraph becomes K  > n and n > K, whereK  ≡ (K + KW )(K + KW + 1)/2. We could now increase KW to the extent
that K  > n condition is satisﬁed9 . In summary, we now have the control
over the number of additional orthogonal variables, KW , so that the above
two conditions are met, while the e vector is unchanged.

8

It is well known that the OLSE of ˛ for the transformed regression, Z  y = Z  X˛ +Z ›, is identical to the OLSE for the original equation. Incidentally, Amemiya op.cit.
showed that the GLS estimator for the transformed model is more eﬃcient than the
OLSE.
9For instance, take the previous numerical case of n = 50 and K = 3. If we set
7 < KW < 47 then the conditions are satisﬁed

10

3.3.3

HCCM Information Accomodated in Candidate Density

We have thus shown that using Z instead of X as the regressor matrix inequation (1), the ω vector is always identiﬁable. Let the OLSE of σ 2 ω beσ 2 ω̂. “σ 2 ” has no hat on it since in our Bayesian MCMC, σ 2 is alreadygenerated in the previous step by the Gibbs sampler scheme. We turn to an
M-H algorithm, to generate ωi . For this purpose, we use a truncated normal
as our candidate generating density:ωi ∼ N (ω̂i , c̃ × ω̂i2 )1ωi ∈(.25, 12.5) ,(13)

for i = 1, . . . , n where c̃ is a constant that has been determined in thefollowing way. In the footnote that discusses inverse gamma distribution(footnote 7), we see that Var(x) = E(x)2 (α − 2)−1 for x ∼ IG(α, β). Hence,
if we were to let α a little over “2” so that the variance would not vanish, thenalmost any value of c̃ suﬃces. The value ought to be such that the average
acceptance rate in the independence chain M-H to be 20% to 50%. The value
we came up with is c̃ = 15 after some experiments, and the acceptance rate
for such candidate generating density function, in our experiment, is 33%,
which is reasonable10 . For r ≥ 2 sequence, equation (13) has been emplyed
as the candidate generating density.

3.4

Hybrid of Gibbs Sampler and M-H

Let us describe our MCMC sequence in its entirety. For the parameter vectorθ = (β  σ ω  ) , we set the initial value θ (0) = (β̂ s2 ιn ) to generate the ﬁrst
round θ (1) . Note here that we ﬁrst draw β (1) , hence β̂ is not really needed inθ (0) . What we actually need is θ(0)
. The sequence of our MCMC continues−β10In M-H, an acceptance rate between 20% and 50% is regarded as reasonable, e.g., see
Koop (2003, p.98). In fact in a series of experiments, we found that the average acceptance
rate of 21% when c̃ = 10, 49% when c̃ = 25, and 58% when c̃ = 35. We may add that
despite these diﬀerences in the average acceptance rates, however, the results are fairly
robust.

11

in the way speciﬁed by equations (6), (8) and (13) for r ≥ 2, combining theGibbs sampler and M-H. The details of monitoring and diagnostic checking
concerning this MCMC sequence, is given in the next section.

4

Numerical and Empirical Investigation

We ﬁrst generate a data set to make sure that our Bayesian procedure works.
After that, we move on to use a Japanese stock price data set to showthat our Bayesian method can be used to estimate stock return volatilities
without assuming any structure on them a priori.

4.1

Numerical Illustration

For the generation purpose, we let the parameters of our model be β =(1, 3, −2) ; σ 2 = 10; sample size n = 50; X = (ιn , x2 , x3 ) ∼ n × 3, whereιn is an n × 1 vector of all one’s, x2 ∼ N (3, In ), and x3 ∼ N (6, 22 In ). Asto ω, we need to specify total of 50 elements. We generated ωi ’s using an
uniform random number generator. Assigned ωi ’s will be given in a chartthat compares them to the Bayesian posterior mean values, later (see Figure
1 below). Using these prespeciﬁed set of parameter values, we generated
data on y. This is used to generate W matrix to form Z = (X W )11 . By
regressing y on Z, we obatin one set of e ∼ n × 1, and hence a ĥ ∼ K  × 1
vector. We now have one sample of y, X and ĥ.
In the MCMC computation, we drew total of 30,000 θ’s. After someexperiments, using such diagnostics devices as Geweke’s (1992), we decided
to throw away the ﬁrst 10,000 simulated values as burn-in’s, and took everyﬁfth drawn values after that to lessen the possible serial correlation inherent
in MCMC. The average acceptance rate in the M-H step was 33%. We
11In our GAUSS programming, we used the “NULL(·)” command to obtain W ∼n × KW that is orthogonal to X. NULL command uses the QR decomposition to obtain
the orthogonal complement of X.

12

Figure 1: Posterior Means of the 50 ωi ’s

now have a set of 20, 000/5 = 4, 000 simulated values for each of σ 2 andωi ’s (i = 1, . . . , 50) available. Since giving 50 posterior pdf’s would notcontribute to the understanding of our procedure, we decided to present
posterior means of the ﬁfty ωi ’s, and this is given in Figure 1.
Fifty assigned ωi values are given using a solid line, and posterior means
using a dotted line. We see that there are some discrepancies between thetrue (or given) value and the posterior mean value, however, mostly they
are close. Their closeness may be observed more clearly for the ﬁrst ﬁve ωi ’sof Figure 1, in Figure 2. In particular, posterior means of ω1 to ω3 , i.e.,E(ω1 |data) to E(ω3 |data), are seen to follow randomly assigned values very
closely. To further this argument, we used some pseudo R–squared measurese.g., given in Wooldridge (2002, pp.465-466), as measures of closeness to the2assigned (or true) ω. Let one of such pseudo R–squared be R̃[E(ω |data)] .
2

In all cases that we computed, we found R̃E(ω |data) to be signiﬁcantly22larger than R̃homogeneous , where R̃homogeneous is the pseudo R–squared13

Figure 2: Posterior Means of the First Five ωi ’smeasure for the homogeneous case, i.e., ω = “a constant vector.” This shows
our E(ωi |data)’s closeness to the assigned ωi .
We depicted the marginal posterior pdf’s of the ﬁrst two σ 2 ωi ’s in Figure
3. We see that the two marginal posterior pdf’s seem to be equipped withtypical characteristics of a variance pdf, i.e., a gamma like shape with ahump close to the origin and skewed to the right. We emphasize that theshapes have been obtained without any prior shape information on them.
This fact may indicate an appropriateness of our MCMC simulator.
We next examine the posterior standard deviations computed from themarginal densities of ωi ’s. The 50 posterior standard deviations are given
in Figure 4. We observe that they tend to ﬂuctuate a lot. For instance,
posterior standard deviation of ω1 and ω2 are slightly less than “1.8” whilethat of ω4 is close to zero. An interesting observation may be made, when
we compare Figure 4 to Figure 1 or to Figure 2. There is a tendency that14

Figure 3: Marginal Posterior’s of σ 2 ω1 and σ 2 ω2Figure 4: Posterior Standard Deviations of the 50 ωi ’s

15

posterior means and standard deviations, move in the same direction. Thismay be interpreted to come from the fact that marginal posterior pdf of ωi ’s
are in gamma like shape12 .

4.2

Empirical Investigation: Japanese Stock Returns

4.2.1

Theoretical ModelIn this section, we use a panel data set to show that our procedure works
well in practice. Saying “panel data set” may be far stretched, however,since we do not implement, for example, random eﬀects model. We carryout a two-step time series to cross section regressions, in a way similar to
the Fama-Macbeth procedure (see e.g., Cochrane (2001 p.244)). We useddata on the daily stock prices of ﬁfty Japanese pharmaceutical/ biomedical
venture capital companies. To obtain an excess return series, we used 10year Japanese Government Bond (JGB) rate for the risk free rate. For themarket return, we used TOPIX. The sample period is from May 6, 2005 toApril 28, 2006, hence the sample size is 245 in total. We obtained stock
return data and the JGB data from Yahoo Finance and Nikko Financial
Intelligence web site, respectively.
We begin with a multifactor multivariate return generating equation:


R = (ιT F )

α
B

+ ε,

(14)

R = (R1 · · · RN ) ∼ T × N is a T period excess returns for N ﬁrms,
α = (α1 , . . . , αN ) ∼ 1 × N vector of constants,
B = (β 1 · · · β N ) ∼ F × N is a matrix of beta’s,
F = (f 1 · · · f F ) ∼ T × F is a matrix of F factors,
ε = (ε1 · · · εN ) ∼ T × N matrix of error terms,
N is the number of stocks, and T is the time series sample size. When F = 1
where

the above reduces to a single factor return generating equation:


R = (ιT f )
12

α
β

+ ε,

(15)

Suppose a gamma distributed random variable γ has E(γ) = θ, then its standard√
deviation becomes SD(γ) = θ/ η, where η is the scale parameter.

16whereF =f
B=β

is a T × 1 vector of one factor,
is a 1 × N vector of beta’s.

Equation (15) is just a set of N time series regressions. We obtain an OLSE
of β , β̂ ∼ 1 × N from equation (15). Deﬁne sample mean of R to be an N
dimensional vector R̄, we then obtain a cross sectional regression modelR̄ = β̂ λ + u,

(16)

1 
R ιT = (R̄1 · · · R̄N ) ∼ N × 1 vector of average excess returns,T
λ ∼ 1 × 1 scalar is a risk premium associated with the factor f ,u ∼ N × 1 is a vector of pricing errors.
Equation (16) is the one factor type CAPM without an intercept term given
where

R̄ =

in Cochrane op.cit., p.235, among others13 .
In this section so far, we have used a set of notations that are common in
empirical ﬁnance, and in this paper it corresponds to equation (1). We, thus,
need to clarify the notational correspondences between the ones used so far
in the current section, and in the previous subsections. The correspondences
are given in below.

R̄ ∼ N × 1β̂ ∼ N × 1
u∼N ×1λ ∼ scalar
F = 1 is the number of factors
N the number of stocks

4.2.2(corresponds to y in (1),(corresponds to X in (1)),(corresponds to ε in (1)),(corresponds to β in (1)),(corresponds to K in (1)),(corresponds to n in (1)).

Bayesian Test of Hetoroscedasticity

Since we do not know whether our data set exhibits a heteroscedasticity ornot, we carried out a Bayesian test of heteroscedasticity a la Lancaster(1994,
p.149). In equation (1) or equivalently in the present case equation (16),13This speciﬁcation is found e.g., in Cochrane op.cit. equation 12.10. We have regressedwith an intercept term, and the OLSE of it is 0.003 (0.048) and the slope estimate is 0.104(0.062), where the ﬁgures inside the parentheses are estimated standard errors. Without
an intercept term, the slope estimate is 0.108 (0.024), and there is very little diﬀerence
whether we include an intercept term or not.

17

under the null hypothesis of homogeneity, H0 , we have ε ∼ N(0, σ 2 I N ).
Hence, letting εi = R̄i − β̂i λ, we have E(ε2i − σ 2 )2 = E(ε4i − σ 2 ε2i + σ 4 ) = 2σ 4 ,under H0 , since our normality assumption on εi calls for E(ε4i ) = 3σ 4 . Lancaster op.cit., thus, proposes T (λ, σ 2 ) =N
22 2
i=1 (ε̂i − σ )N− 2σ 4 , whereε̂i = R̄i − β̂i λ̂ and λ̂ is an MCMC generated value, to be used in a Highest
Posterior Density Interval (HPDI) test. This test would reject H0 if the massof the T statistic distribution does not contain “zero,” i.e., homoscedasticity14 .
It is easy to draw λ and σ from their respective marginal posterior pdf,
since each pdf boils down to an univariate normal and an IG, respectively.
We will not discuss the details of the MCMC procedure involved in simulating T under H0 . Simulated T statistic is shown in Fig. 4. Our calculations show that for the distribution of T , P (T < 0) = 0.016, which issmaller than the usual signiﬁcance level even if we multiply it by 2, i.e.,
0.016 × 2 = 3.2% < 5%, say. We would safely reject H0 . We also carried out
several well-known sampling theory heteroscedasticity tests. They all reject
H0 .
4.2.3

Empirical Estimation of ωi ’s

The MCMC simulation procedure, here, is exactly the same as the one
discussed in the numerical illustration of section 4.1. Posterior means andstandard deviations are given, respectively, in Figures 5 and 6 in below.
Since we have conﬁrmed that our method works in section 4.1, we need to14
Lancaster op. cit. does not call T , a test statistic, however. He rightfully regards T
as a “checking” device. For one thing, the empirical distribution obtained from simulated
T is not the marginal posterior pdf under H0 . In other words, we (or Lancaster op cit.
too) did not compute the posterior pdf of the T statistic. Secondly, in a sampling theory
framework, obviously plimT = 0 under H0 as N → ∞. This implies that whether T > 0
or T < 0 is not known under the alternative hypothesis as N → ∞. It seems like we ought
to take a 100(1 − α)% HPDI of the T statistic, and reject H0 if it does not contain zero,
where α is a signiﬁcance level. This is in essence, what we did.

18

Figure 5: Marginal Posterior pdf of the T Statistic

Figure 6: Marginal Posterior Means of the σ 2 ωi ’s

19

Figure 7: Marginal Posterior Standard Deviations of the σ 2 ωi ’s

make sure that our results may be reasonable as an empirical investigation.

In the U.S. and world wide, starting with a seminal paper by Banz
(1981) many observed “size eﬀect” that is the smaller the company is thehigher its return15 . This phenomenon could be naively associated to themean-variance eﬃciency to yield a thesis that says smaller companies are
expected to be more riskier, i.e., the smaller the size, the larger the mean
and volatility of returns. In this section, we intend to investigate whether
larger (smaller) size companies have smaller (larger) ωi ’s.
For the ﬁfty stock data, we now collected the market value of its equity(MVE) data from Yahoo Finance (in one billion yen)16 . We then drew agraph with the MVE on the horizontal axis, and the σ 2 ωi ’s on the vertical
15

For the size eﬀect in Japan, see e.g., Chan and Chen (1991), among others.
Berk (1997) among others, however, examined the so-called size eﬀect and proposedthat a size of a company should not be measured by MVE but some other variables suchas sales.
16

20

Figure 8: Volatility versus Size: σ 2 ωi ’s and MVE

axis. This is shown in Figure 7. The downward sloping solid line in Figure
7, is the OLS estimated line. This ﬁgure clearly shows the larger the MVE,
the smaller the volatility as measured by σ 2 ωi . We further selected two
stocks (1) that has large MVE and small σ 2 ωi , and (2) that has small MVE
and large σ 2 ωi , to see what the marginal posterior pdf’s of σ 2 ωi of thesecompanies look like. Speciﬁcally, we chose Taisho Pharmaceutical Co., Ltd.
for (1), and Site Support Institute Co., Ltd. for (2). They are given in
Figures 8 below. Notice that the two charts have diﬀerent horizontal axis
scale. The smaller MVE stock has very large volatility (see the lower chart)
compared to the that of the larger MVE stock (see the upper chart). The
two pdf’s have quite reasonable shapes. We conclude that our Bayesianestimation of volatility supports the view that the smaller the size of thestock, the larger is the volatility.

21

Figure 9: Marginal Posterior pdf of σ 2 ωi for Taisho Phamaceutical Co. andSite Support Institute5

Concluding Remarks

In this paper, we proposed a fully Bayesian method to estimate regression
error term scedastic function ω without assuming any structure on it. One
unique feature of our methodology is in using the Eicker-White result ina candidate generating function in our M-H step of MCMC. We have supplied a numerical example, and also an empirical example to show that our
method works reasonably. Finally, as demonstrated in our empirical research, obtaining W such that W  (y X) = 0 is no problem since W is notfound as an empirical data but rather computer generated data.
One possible arbitrariness might remain in this paper. That is, a reader
might feel uneasy about our setting bounds for ω prior pdf. As we suggested
in footnote 8, we may adopt a hierarchical Bayes model.
Acknowledgements

22This research is supported in part by the Ministry of Education and Culture,
Grant in Aid #19530178. We thank Hiroki Tsurumi of Rutgers University,
and Takeshi Amemiya of Stanford University for incisive comments andsuggestions for earlier drafts of this paper. We are responsible for all errorsthat inadvertently remain in this paper, however.

References[1] Amemiya, T., 1983, Partially generalized least squares and two-stage
least squares estimators. Journal of Econometrics, 23, 275–283.
[2] Amemiya, T., 1985, Advanced Econometrics, Harvard University Press.
[3] Banz, R.W., 1981, The relationship between return and market valueof common stocks. Journal of Financial Economics, 9, 3–18.
[4] Belsley, D.A., 2002, An investigation of an unbiased correction for heteroskedasticity and the eﬀects of misspecifying the skedastic function.
Journal of Economic dynamics and Control, 26, 1379–1396.
[5] Berk, J.B., 1997, Does size really matter? Financial Analysts Journal,
53, no.5, 12–18.
[6] Chan, K.C. and Nai-Fu Chen, 1991, Structural and return characteristics of small and large ﬁrms. Journal of Finance, 46, 1467–84.
[7] Chang, E.C. and S. Dong, 2006, Idiosyncratic volatility, fundamentals,
and institutional herding: Evidence from the Japanese stock market.
Pacific-Basin Finance Journal, 14, 135–154.
[8] Cochrane, J.H., 2001, Asset Pricing, Princeton University Press.

23[9] Eicker, F., 1963, Asymptotic normality and consistency of the leastsquares estimators for families of linear regressions. Annals of the Mathematical Statistics, 34, 447–456.
[10] Geweke, J., 1992, Evaluating the accuracy of sampling-based approaches to the calculations of posterior moments. in Bernardo, J.,
Berger, J., Dawid A., and Smith, A. (eds.), Basyesian Statistics, Volume 4, pp.641–649, Clarendon Press.
[11] Godfrey, L.G., 2006, Tests for regression models with heteroskedasticity
of unknown form. Computational Statistics & Data Analysis, 50, 2715–
2733.
[12] Greene, W.H., 2008, Econometric Analysis, 6th ed., Prentice-Hall.
[13] Koop, G., 2003, Bayesian Econometrics, John Wiley & Sons.
[14] Long, J.S. and L.H. Ervin, 2000, Using heteroscedasticity consistent
standard errors in the linear regression model. American Statistician,54, 3, 217–224.
[15] Lütkepohl, H., 1996, Handbook of Matrices, John Wiley & Sons.
[16] Robinson, P.M., 1987, Asymptotically eﬃcient estimation in the presence of heteroskedasticity of unknown form. Econometrica, 55, 4, 875–891.
[17] Robert, C.P. and G. Casella, 2004, Monte Carlo Statistical Methods,
2nd ed., Springer.
[18] White, H., 1980, A heteroscedasticity-consistent covariance matrix estimator and a direct test for heteroscedasticity. Econometrica, 48, 817–
838.
24

[19] Wooldridge, J.M., 2002, Econometric Analysis of Cross Section andPanel Data, MIT Press.

Appendix : Derivation of equation (11)
H is deﬁned by H = X 0 ΩX. Our purpose in this appendix is to derive
equation (9) in the text, ı.e., h ≡ vech(H) = Xn ω. To this eﬀect, we notethe following series of lemmas.
Lemma 1 . vec(H ) = (X  ⊗ X  )vec(Ω) ∼ K 2 × 1.
proof : Just use the usual vec(ABC) = (C  ⊗ A)vec(B) formula. See,e.g., Lütkepohl (1996, p.97).
Lemma 2 . X  ⊗X  = [(x1 ⊗x1 ) . . . (x1 ⊗xn )] . . . [(xn ⊗x1 ) . . . (xn ⊗xn )] ∼K 2 × n2 .
proof : Note X  = (x1 . . . xn ) ∼ K × n, and the result follows.
Lemma 3 . vec(H) = [(x1 ⊗ x1 ) . . . (xn ⊗ xn )]ω, where ω = (ω1 . . . ωn ) ∼n × 1.
proof : Note.
.
.
vec(Ω) = vec(diag(ω)) = (ω1 0 · · · 0 .. 0 ω2 0 · · · 0 .. · · · .. · · · ωn ) ∼ n2 ×1,
where ω  = (ω1 , . . . , ωn ), “diag(ω)” denotes a diagonal matrix with thevector ω in its diagonal. Note that vec(Ω) is a vector of size n2 × 1,
and picks up the columns (xt ⊗ xt ) only for t = 1, . . . , n from X  ⊗ X  .
Lemmas 1 and 2, in conjunction with the above proves the lemma.
Lemma 4 . xt ⊗ xt =vec(xt xt ) ∼ K 2 × 1.

25

proof : Use the vec(AB) = (B  ⊗ A)vec(I p ) relationship, where p is thecolumn dimension of A (see e.g., Lütkepohl (1996, p.97)). In thepresent case p = 1, hence, vec(I p ) is equal to a scalar one.
We now deﬁne a matrix that converts a vec(A) type vector to a vech(A)vector.
Deﬁnition of LK  matrix : Let K  ≡ K + (K − 1) + . . . + 1 = 12 K(K + 1)
and deﬁne
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
LK  = ⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣

IK

0

0
.. I k−1.

0

0 0.. .. I k−2. .
··· ···
0 ··· 0 1

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥ ∼ K  × K 2.
⎥
⎥
⎥
⎥
⎥
⎥
⎦

In essence, the diagonal submatrices of LK  shrink from I K , I K−1 . . . I K−j . . .
to I K−(K−1) = 1 but each submatrix is preceded by a zero matrix 0j ,where 0j ∼ (K − j) × j. Such LK  converts vec(A) to vech(A) for a
K × K square matrix A as follows.
Lemma 5 . Let A ∼ K × K square, then vech(A) = LK  vec(A) ∼ K  × 1.
proof : Let A = (a1 . . . aK ) ∼ K × K and aj ∼ K × 1. Also let, forinstance,aj(−1,2) = (a3j a4j · · · aKj ) ∼ (K − 2) × 1,

26

i.e., aj(−1,2) deletes the ﬁrst two elements from aj . Then,⎛⎞a1⎜ ...... ⎟
⎜
⎟
⎜
⎟
0
⎜
⎟
⎜ a2(−1) ⎟⎛⎞⎜ ...... ⎟
⎜
⎟
a1
⎜
⎟
⎜ . ⎟
0
⎜
⎟
⎟=L ⎜.
LK  vec(A) = LK  ⎜⎟ = vech(A).
⎝ . ⎠
K ⎜
⎟
0
⎜
⎟
aK
⎜ a
⎟
⎜ 3(−1,2) ⎟⎜ ...... ⎟
⎜
⎟
..
⎜
⎟
⎝
⎠
.
aKK
Lemma 6 . vech(xt xt ) = LK  vec(xt xt ).
proof : Omitted.
We are now, in a position to derive the vech(H) = Xn ω relationship.
Use Lemmas 2 to 4 on the Lemma 1 relationship to obtainvec(H) = [vec(x1 x1 ) . . . vec(xn xn )]ω.
Further, use Lemmas 5 and 6 to obtain the following.
vech(H) = LK  vec(H)ω
= [LK  vec(x1 x1 ) . . . LK  vec(xn xn )]ω
= X n ω.
where Xn = [vech(x1 x1 ) . . .vech(xn xn )] ∼ K  × n.

27