Discussion Paper No. 341
DIRICHLET PRIOR FOR ESTIMATINGUNKNOWN REGRESSION ERRORHETEROSKEDASTICITY
Hiroaki Chigira
Tsunemasa Shiba
December 22, 2015

TOHOKU ECONOMICS RESEARCH GROUPGRADUATE SCHOOL OF ECONOMICS AND
MANAGEMENT TOHOKU UNIVERSITY

27-1 KAWAUCHI, AOBA-KU, SENDAI,
980-8576 JAPANDIRICHLET PRIOR FOR ESTIMATING UNKNOWN
REGRESSION ERROR HETEROSKEDASTICITY∗
Hiroaki Chigira† and Tsunemasa Shiba‡
12/2015

Abstract
We propose a Bayesian procedure to estimate heteroskedastic variances of the regression error
term ω, when the form of heteroskedasticity is unknown. The prior information on ω is basedon a Dirichlet distribution, and in the Markov Chain Monte Carlo sampling, its proposal density
parameters’ information is elicited from the well-known Eicker–White Heteroskedasticity Consistent
Variance–Covariance Matrix Estimator. We present an emprical example to show that our scheme
works.

key words
Dirichlet prior, Eicker–White HCCM, informative prior pdf’s, MCMC,JEL Classification
C11, C13∗

We thank Hiroki Tsurumi, Takeshi Amemiya, and Toshiaki Watanabe for incisive comments and suggestions forearlier drafts of this paper. We are responsible for all errors that inadvertently remain in this paper, however.
†Department of Economics, Tohoku University, 27-1, Kawauchi, Aoba-ku, Sendai-shi 980-8576 Japan. e-mail: hchigira@econ.tohoku.ac.jp‡
Graduate School of Finance, Accounting & Law, Waseda University, 1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027,
e-mail: tshiba@waseda.jp, and Department of Economics, Hitotsubashi University.

–1–Dirichlet Prior for Estimating Unknown Regression Error Heteroskedasticity
Hiroaki Chigira and Tsunemasa Shiba

1

Introduction

We propose a Bayesian approach to estimate heteroskedastic parameters of regression error variancesthat are of unknown form, using Dirichlet prior pdf. As Amemiya (1985, p.199) points out, thecrucial ω 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 ω. Eicker (1963) andWhite (1980) independently developed a well-known consistent variance-covariance matrix estimator(“HCCM” hereafter) for the OLS regression coeﬃcient estimator. We use HCCM information toformulate proposal density of a Metropolis-Hastings (“M-H” hereafter) algorithm in Markov Chain
Monte Carlo simulation. Unidentiﬁability of ω poses no problem to us. As in Amemiya (1985, p.199)
we use an orthogonal regression that circumvents possible underidentiﬁability of ω, this problem. Also,a Dirichlet prior on ω should make it identiﬁable in a Bayesian context.
The trend in the HCCM literature seems to be in the ways to improve ﬁnite sample performanceof tests of the linear restriction(s) on the coeﬃcient vector, e.g., Long and Ervin(2000) and Godfrey(2006), among others. Our focus in this paper is in the direct estimation of the elements of ω. There
are papers that deal with statistical inferences of regression coeﬃcients, when the skedastic functionof the error term is unconstrained. Robinson (1987) assumes it to be a function of regressors, andderives an GLS estimator that is more eﬃcient than the existing ones. Our Bayesian estimation of ω
will help to sharpen posterior density of β and/or lead to a better predictive density. It may also leadto more eﬃcient estimator of β in terms of asymptotic theory framework as well.
In Bayesian econometrics, starting with a seminal work by Geweke (1993) there is a homoskedastic
Student-t regression model derived from normal linear regression (“NLR” hereafter) with a particular
set of Gamma priors for heteroskedasticity parameters. This model has been introduced in such books
as Koop (2003), Geweke(2005) and Greenberg (2013), among others, and is now a popular Bayesianmodel. We will compare our model to this model using the Deviance Information Criterion (“DIC”
hereafter).
We shall discuss estimation of ω, where ωi is the ith volatility. In order to access an option
pricing, we ﬁrst need to come up with a reasonable estimate of volatility. Our estimation of ω needs
1

The ω vector has in its elements, all the normalized diagonal elements of variance-covariance matrix of the regression
error term. The normalization rule for the matrix is given just below equation (4).

–2–

no parametric model for volatility process such as the GARCH model, since we only use HCCM
information. If we wish to estimate a volatility 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 atheoretical background. Moreover, when it comes to cross section data, historical volatility calculation
breaks down for obvious reasons.
After assuming a usual prior density for the parameters in the regression model, we may obtain ajoint posterior density. The usual parameters, e.g., (β, σ 2 ), may be easily simulated using the Gibbs
sampling. It is in the simulation of the elements of ω that we use the HCCM. We use results from
HCCM to form the proposal 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 a Gibbs sampler. We propose to use the Eicker–White result to simulate ω bya M-H scheme. In sections 4 numerical illustrations to compare our methodology to homoskedastic t
distiruted error term model, are given. We us DIC to this eﬀect.

2

Model and Assumptions

2.1

Likelihood

Let an NLR model with heteroskedastic error term be
yi = x′i β + ui , (i = 1, . . . , n)(1)where yi ∼ dependent variable, xi ∼ K ×1 non stochastic explanatory variables, β ∼ K ×1 coeﬃcients,
and the properties of regression error term u be
ui | ωi , σ 2 ∼ N (0, σ 2 ωi ) .

(2)

Our single likelihood function for yi has the following normal distribution
yi | xi , β, ωi , σ 2 ∼ N (x′i β, σ 2 ωi ) .
We may use the following two notations: ωi =(3)

1
, where λi is the precision. Geweke (1993) uses ωi ,
λiwhereas Geweke (2005) and Greenberg (2013) use λi . (1) in an matrix form becomesy = Xβ + u ,

–3–

(4)

where y ∼ n × 1 of yi ’s, X ∼ n × K matrix of stacked up x′i , u ∼ Nn (0, σ 2 Ω), and Ω = diag(ω) =diag(ω1 , . . . , ωn ) ∼ n × n withn∑

ωi = tr(Ω) = n. Note that tr(Ω) = n restriction is often employed in

i=1heteroskedasticity literature, e.g., Greene (2012, p 308). The likelihood function for the whole sample
becomesℓ(y|X, θ) ∝ σ −n (

n
∏

(

−1

ωi 2 ) exp

i=1−1

where we noted |σ 2 Ω| 2

= σ −nn∏)−1(y − Xβ)′ Ω−1 (y − Xβ)
2σ 2

,

(5)

−1

ωi 2 . On completing squares on β above, we obtain

i=1ℓ(y|X, θ) ∝ σ −n (

n
∏

i=1

(

−1

ωi 2 ) exp)−1′(νs2 + (β − β̂)′ X̃ X̃(β − β̂) ,
2
2σ

(6)

where ν = n − K, νs2 ∼ the sum of squared residuals from the regression of ỹ on X̃, ỹ = Ω−1/2 y,′′X̃ = Ω−1/2 X, and β̂ = (X̃ X̃)−1 X̃ ỹ = (X ′ Ω−1 X)−1 X ′ Ω−1 y is the GLS estimator of β. (6) turns
out to be useful for simulating β since it is in a multivariate normal form in β.

2.2

Prior and Posterior

Let θ = (β ′ , σ 2 , ω ′ )′ ∼ (K + 1 + n) × 1 be our parameter vector. We assume prior distribution on θ asπ(θ) ∝ π(β) π(σ 2 ) π(ω) ,(7)

where the three arguments of θ are independent. Individual prior for β and σ 2 areβ ∼ NK (β 0 , B 0 ), σ 2 ∼ IG(

α0 δ 0
, ),
2 2

(8)

where β 0 ∼ K × 1, B 0 ∼ K × K, α0 and δ0 are hyper parameters in the prior pdf’s that are assumed
to be known, and IG(·) denotes an inverted gamma distribution.
We note that ω is an n dimensional continuous multivariate random variable vector that satisﬁes
the tr(Ω) =∑i ωi = n restriction. The best suited prior to this regard, is obviously Dirichlet with its

hyper parameter values the same for all i = 1, . . . , n. This way, we may eﬀectively represent unknown
heteroskedasticity structure and satisfy tr(Ω) = n restriction. If we employed Gamma prior, e.g.,
Geweke (1993, 2005) and Greenberg (2008), among others, then we are in eﬀect imposing a certain
structure in the heteroskedasticity.
Since a Dirichlet has the property that its elements add up to one, not n, we cannot place a
Dirichlet prior directly on ω. Instead, we shall assumeω̃ ∼ D(η 0 ),
–4–

(9)

where ω̃ =1
ω, and D(η 0 ) denotes a Dirichlet distribution with a parameter vector η 0 = (η01 , . . . , η0n )′ ∼nn × 1. The values of η 0 will be given later in this paper. The assumption on ω̃ is tr(Ω̃) = 1, thustr(Ω) = n. If we make a transformation from ω̃ vector to ω vector, we arrive at our prior distribution
on ω that resembles to Dirichlet aside from normalizing constant:
(
∑

)

(

)

∑

Γ
η0i
Γ
η0i n
)
n (
n−1−η̃0∏∏ωh η0h −11iiπ(ω) =
ωℓ η0ℓ −1 ,
= n
n
n
∏
n ∏
n
Γ(η0j ) h=1
Γ(η0j ) ℓ=1
j=1

where η̃0 =

n
∑

(10)

j=1

η0h . Hence, its kernel is given by

h=1

n
∏

π(ω) ∝

η −1

ωi 0i

.

(11)

i=1

Given (8) and (11), our joint prior for θ becomes(−1qπ(θ) ∝ exp
2 β

)(

1
σ2

) α0 +1
2

(

−δ0
exp
2σ 2

)∏
n

η

ω̃h 0h

−1

,

(12)

h=1

where qβ = (β − β 0 )′ B −10 (β − β 0 ). Finally, the posterior π(θ|y, X), is obtained by combining (6) and(12) as π(θ|y, X) ∝ ℓ(y|X, θ)π(θ) to obtainπ(θ|y, X) ∝ σ −n (

n
∏

−1

ωi 2 ) exp

i=1

(

(
)
(
) ( ) α0) n−( 2 +1)
−ψ
−δ0 ∏
−1
η −1
2
q
exp
exp
σ
ω̃h 0h ,β222σ22σ h=1(13)′

where ψ = νs2 + (β − β̂)′ X̃ X̃(β − β̂). Note that in θ, ψ depends on β and ω, while qβ depends on β.

2.3

Student-t homoskedasticity model and our model

The main purpose of our model is in Bayesian estimation of the n elements in ω without assuming
a prior structure for it. On the other hand, Geweke’s (1993) Student-t homoskedasticity model, also
given in Geweke (2005) and Greenberg (2012), among others, is primarily concerned with estimating
β coeﬃcient vector in heteroskedastic NLR model.
Suppose a heteroskedastic regression disturbance ui | λi , σ 2 ∼ N (0, σ 2 /λi ). If we assume a Gamma
with an identical shape and scale parameters, ν = ν1 = ν2 on a prior for λi , then the resultant ui | λi , σ 2
distribution becomes a fat-tailed homoskedastic Student-t with the parameters (ν, 0, σ 2 ).
If there is a compelling need for such Gamma prior on λi , then the above equivalence of (1) NLR
with heteroskedasticity, and (2) homoskedastic Student-t should be of a great value. What would
the eﬀect be if there is only one parameter value ν1 = ν2 = ν in the Gamma distribution? Such
–5–

Gamma random variable would have E(x) = 1 always, and its pdf becomes concentrated abround it
as ν gets large. We may thus conclude that the single parameter Gamma distribution assumption,is rather peculiar. In view of the above conclusion, we suggest that we depart from λi ∼ Gamma
distribution assumption, and adopt more reasonable prior for ui heteroskedasticity. Notice, however,that even if a Gamma prior G(ν1 ν2,), on λi is assumed, then n λi ’s are, again, generated from one
2 2

single Gamma. This is in eﬀect assuming a particular structure on the heteroskedasticity of ui ’s. If
the interest of the Bayesian analysis, is in ﬁnding the parameter of the structure, ν = ν1 = ν2 , thenthe Gamma distribution assumption may be justiﬁed. But if we are interested in estimating each
σi2 = σ 2 /λi = σ 2 ωi then we need something other than a scalar Gamma assumption. As we developedin the previous subsection, we employ a Dirichlet prior to this end.
Dirichlet prior on λi or ωi is suitable for Bayesian estimation of heteroskedastic variance parametersthat have unknown structure, on two grounds. First, it is free of restrictions from the small number ofparameters that govern the entire shape of the prior pdf. In Geweke’s Gamma pdf prior for λ ∼ n × 1,
ν is the only parameter of the distribution. If the prior pdf of the λ vector were to be of unknown
structure, it should have n parameters. Secondly, Dirichlet distributed random variables, xi ’s arecontinuous and satisfy the
setup, where

3

∑

n
∑

xi = 1 constraint by construction. This is a welcome constraint to ouri=1

ωi = n restriction, needs to be satisﬁed a priori.

MCMC Simulation of θ

We use notations “θ −ϑ ” to denote θ less ϑ hereafter. For instance, “θ −β ” implies θ −β = (σ 2 , ω ′ )′ ∼(1 + n) × 1.

3.1

Gibbs Sampler for β and σ 2

As shown in the two remarks below, tractable fully conditional posteriors of β and σ 2 may be obtained.
On the other hand we need to implement an Independence Chain M-H algorithm for simulating ω.
Remark 1.

Fully conditional posterior of β is given byβ| θ −β , y, X ∼ NK (β 1 , σ −2 B 1 ) ,
′

(14)

′

where B 1 = (X̃ X̃ + (σ −2 B 0 )−1 )−1 , β 1 = B 1 φ, and φ = X̃ X̃ β̂ + (σ −2 B 0 )−1 β 0 .
Proof From the joint posterior (13), conditional posterior of β becomes
(

)

(

−1
−ψ
exp
qβ
π(β| θ −β , y, X) ∝ exp(
2
2σ
2
–6–

)

∝ exp(−1Aβ ) ,
2σ 2′where Aβ = (β − β̂)′ X̃ X̃(β − β̂) + (β − β 0 )′ (σ −2 B 0 )−1 (β − β 0 ). On completing squares for β, Aβ
becomes′′−1
Aβ = (β − β 1 )′ B −1
+ σ −2 B 0 ]−1 (β̂ − β 0 ) .
1 (β − β 1 ) + (β̂ − β 0 ) [(X̃ X̃)Hence,()

1π(β| θ −β , y, X) ∝ exp − (β − β 1 )′ (σ −2 B 1 )−1 (β − β 1 )2(15)The right hand side of (15) may be used to simulate β from a multivariate normal.
Remark 2.

Fully conditional posterior of σ 2 is in Inverted Gamma:(σ | θ −σ2 , y, X ∼ IG
2n + α0 ψ + δ 0
,
2
2

)

(16)

Proof From the joint posterior (13), conditional posterior of σ 2 becomesπ(σ 2 | θ−σ 2 , y, X)∝ σ −n exp( 2 )−)(α0ψ ( 2 )−( 2 +1)
−δ0
− 2 σ
exp
2σ
2σ 2

( n+α0
2

∝ σ

3.2

(

)

+1

(

(ψ + δ0 )
exp −
2σ 2

)

)

.

Independence Chain for ω

We now turn to ω simulation. From the joint posterior, (13), we haveπ(ω| θ −ω , y, X) ∝ (n∏(η ∗ −1)
ωi 0i ) exp(

i=1∗ = η −
where η0i
0i

−ψ
2σ 2

)

,

(17)

1
> 0 for i = 1, .., n in order for η0∗ vector to make sense as a Dirichlet parameter.
2

Let
Aω =

n
∏
i=1

(η ∗ −1)
ωi 0i

(

−ψ
and Bω = exp
2σ 2)

hence π(ω| θ −ω , y, X) = Aω Bω . Obviously Aω is a kernel of D(η0∗ ). On the other hand Bω certainly
looks like a Nn (Xβ, σ 2 Ω), however, as a kernel of ω, Bω is not of any known form.
We shall use an Independence Chain M-H simulator for ω. Since Bω is not going to give any clue
for a proposal density, we use information contained in Aω ∼ Dirichlet distribution, for our proposal
density. Particular value of the parameter vector in the proposal density, will be discussed later. In
the following, we shall ﬁrst give an outline of our M-H strategy.

–7–

3.2.1

M-H Acceptance Probability of ωLet ω (r) be the “rth” current value in the chain. Then the acceptance probability of ω (′) given ω (r)would be

(

α(ω

(r)

π(ω (′) , θ −ω | y, X) f (ω (r) ), ω ) = min 1,
π(ω (r) , θ −ω | y, X) f (ω (′) )

)

(′)

,

where f (ω) represents our proposal density for ω, and it is explained in below. We ﬁrst take up theratio of posterior densities in the acceptance probability. Noting that arguments other than ω cancels
out, it becomes
∏ (

(′)

)η∗ −10ii ωiπ(ω (′) , θ −ω | y, X)=∗ −1∏ ( (r) )η0j
π(ω (r) , θ −ω | y, X)ωj(where ψ (′) = (y − Xβ)′ Ω(′)

)−1

j

(

)

(

),

−ψ (′)exp
2σ 2−ψ (r)
exp
2σ 2

(18)

(y − Xβ), Ω(′) = diag(ω (′) ) and this is not to be confused with a(transpose of Ω. Likewise, ψ (r) = (y − Xβ)′ Ω(r))−1(y − Xβ), Ω(r) = diag(ω (r) ). Suppose that we

employ D(η p ) as our proposal2 , then the ratio of proposal densities above becomes∏ ()(r) ηpi −1i ωif (ω (r) )=∏ ( (′) )ηpj −1 .
f (ω (′) )
ω
j

(19)

j

Results of (18) and (19) may be put together to produce


α(ω

3.2.2

(r)



, ω ) = min 
1,

(′)

(

−ψ (′)
exp
2σ 2
exp

)

∏ (

)η −η∗





(
)
(
)η −η∗ 
.
∏
(r)
(′)
j
0j 
−ψ
ω

2σ 2

(r)

i

ωi

j

j

i

0i

(20)

Proposal Density ParameterIn this subsection, we outline how we specify the Dirichlet proposal density that uses White’s HCCM
information. Notice that our prior for ω = n ω̃ resembles to a Dirichlet, where ω̃ ∼ D(η0 ). We let theparameter vector of the Dirichlet proposal density be η p = c η̂, where c is a scalar tunning parameter,
and information on η̂ is to be extracted from a regression to be explained in below. We brieﬂy discuss
how we obtain η̂.
We ﬁrst regress y on X by the OLS to obtain estimated residual vector eols . We then use it to construct White’s HCCM estimator, Ĥ. Let the vector obtained from Ĥ be ĥ, where ĥ = vech(Ĥ) ∼ K ′ ×

2

The details of our proposal density are given in the next subsection.

–8–

1 and K ′ = 12 K(K + 1). As a regressor matrix to ĥ, consider X n = [vech(x1 x′1 ), . . . , vech(xn x′n )] ∼K ′ × n. Regression of ĥ on X n , i.e.,ĥ = X n σ 2 ω + v ,(21)

where v is some error term vector, will yield an estimator of ω, when K ′ > n and n > K. Whenthese conditions are not met, we may always augment X by W ∼ n × KW such that W ′ (y, X) = 0.
This is essentially ﬁnding an orthogonal comlement of the matrix (y, X). Using such software as R
and GAU SS, among others, we can easily ﬁnd W . Actually we only need W ∼ n × KW such that
KW meets the next two conditions: 21 (K + KW )(K + KW + 1) > n and n > (K + KW ). Since W is
orthogonal to both y and X, OLS estimated regression residual from (21), eols , is the same whether
we used X or W augmented (X, W ) as the regressor matrix. In summary, we may always estimate
ω given data, y and X, using (21). Let the resultant estimator of ω from (21) be ω̂ ols .
It turns out, however, that we cannot simply set ω̂ ols equal to η̂. Our preliminary investigations
revealed that there are cases where the variability of ω̂ ols is too big for it to be η̂. A variability
measure of ω̂ ols could be R =max(ω̂ols i ), where ω̂ols i is the ith element of ω̂. We aim to obtain amin(ω̂ols i )moderate R value, while preserving the ranking, the order or the relative magnitude of the elementsof ω̂ ols that HCCM information conveys. To this eﬀect, we introduce an additional scalar tunningparameter d ≥ 0 to be added to ω̂ols i as η̂i = (ω̂ols i + d), where η̂i is the ith element of η̂. Deﬁne the
new ratio to be f (d) =max(ω̂ols i ) + d. This ratio, f (d) is a function of d ≥ 0, and has the propertymin(ω̂ols i ) + dsuch thatR ≥ f (d) ≥ 1 since f ′ (d) < 0, lim f (d) = 1, and lim f (d) = R.
d→∞

d→0

Thus by an appropriate choice of d, we may make the variability of ω̂ ols to be any desirable level. A
suggested value for d is to set f (d) = n, i.e., sample size n dependent. Solving the equation for d, we
obtaind=max(ω̂ols i ) − n min(ω̂ols i )1≈ max(ω̂ols i ) − min(ω̂ols i ),n−1nwhich is a reasonable value.
In summary our Dirichlet proposal density parameter η p = c η̂ is set to be η̂i = ω̂ols i + d with c andd tunning parameters. One proposed value for d is given in the preceding paragraph. If it happensthat R ≤ n, then set d = 0 and η̂ = ω̂ ols would do the job.
–9–

4

Empirical Investigation: Japanese Stock Returns

In this section, using Japanese stock dataset, we present a comparison of our Dirichlet prior method,
Student-t homoskedasticity model (Geweke(1993)’s), and the usual homoskedastic NLR model, in
terms of DIC.

4.1

One Factor Model and Dataset

We carry out 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 used data on the daily stock prices of ﬁfty Japanesepharmaceutical/ biomedical companies. To obtain an excess return series, we used 10-year Japanese
Government Bond (JGB) rate for the risk free rate. For the market return, we used TOPIX. The
sample period is from May 6, 2005 to April 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 an one factor return generating equation:(

R = (ιT f )

α
β

)

+ ε,

(22)

where R= (R1 · · · RN ) ∼ T × N is a T period excess returns for N ﬁrms, ιT ∼ vector of one’s, f isa T × 1 vector of one factor, α=(α1 , . . . , αN ) ∼ 1 × N vector of constants, β = (β 1 · · · β N ) ∼ 1 × N
is a vector of beta’s, ε=(ε1 · · · εN ) ∼ T × N matrix of error terms, N is the number of stocks, and T
is the time series sample size. Equation (22) is just a set of N time series regressions. We obtain an
OLSE of β , β̂ ∼ 1 × N from equation (22). Deﬁne sample mean of R to be an N dimensional vectorR̄, we then obtain a cross sectional regression model′R̄ = β̂ ϕ + u,
where R̄ =

(23)

1 ′
R ιT = (R̄1 · · · R̄N )′ ∼ N × 1 vector of average excess returns, ϕ ∼ 1 × 1 scalar is a riskTpremium associated with the factor f , u ∼ N × 1 is a vector of pricing errors. Equation (23) is theone factor type CAPM without an intercept term given in Cochrane op.cit., p.235, among others3 .

3This speciﬁcation is found e.g., in Cochrane op.cit. equation 12.10. We have regressed with an intercept term,
and the OLSE of it is 0.001 (0.054) and the slope estimate is 0.109 (0.069), where the ﬁgures inside the parentheses
are estimated standard errors. Without an intercept term, the slope estimate is 0.110 (0.025), and there is very little
diﬀerence whether we include an intercept term or not.

–10–

In this section so far, we have used a set of notations that are common in empirical ﬁnance,
and in this paper equation (23) corresponds to equation (1). We need to clarify the notationalcorrespondences between the ones used so far in the current section, and in the previous subsections.
′

The correspondences are: R̄ ∼ N × 1 corresponds to y in (1), β̂ ∼ N × 1 corresponds to X in (1),u ∼ N × 1 corresponds to ε in (1), λ ∼ scalar corresponds to β in (1), the number of factor in (23) isone and it corresponds to K in (1), N the number of stocks corresponds to n in (1).

4.2

The Three Models Compared

The three models we compare are(i) homoskedastic NLR model,(ii) Geweke’s (1993) model, and(iii) our model in this paper.
We designed the common parameters of the three models to be the same. For instance, all are basedon NLR of equation (23), with mutually independent priors:(ϕ ∼ N (β0 , B0 ), σ ∼ IG
2

α0 δ 0
,
2 2

)

(24)

where we set α0 = δ0 = 10−2 , β0 = 0, and B0 = 104 . The posterior pdf is given byπ(θ|y, X) = π(ϕ)π(σ 2 )π(ω)π(y|ω, X),(25)

where π(ϕ) and π(σ 2 ) are given in (24), while π(ω) depends on the model, and π(y|ω, X) is amultivariate normal based likelihood given in (6)4 . In all three MCMC similations, burn-in is set to
be 10,000 while one apart simulated values are taken out of 8,000, i.e. total of 4,000 simulated values
are obtained.
Let us ﬁrst discuss model (i), homoskedastic NLR Model. Homoskedasticity in terms of π(ω) isrepresented by ω = ιn , i.e., all elements to be one. The model is now a homoskedastic NLR, however,
MCMC is needed since our prior for σ 2 is an informative inverted gamma.
We next take up model (ii), Geweke (1993)’s model. This is an NLR model with heteroskedasticity.
As we stated in section 2.3, Geweke (1993) converted this model to Student-t homoskedastic model.
Unlike our inverted gamma σ 2 prior in (24), he used noninformative σ 2 prior, however. The same4

β in (6) should be replaced by ϕ.

–11–

inverted gamma σ 2 prior as in (24), is used in Koop (section 6.4, 2003). Hence, we decided to follow
Koop (op. cit.)’s MCMC sampler for model (ii). Joint posterior of Geweke (1993)’s model is still (25)
except for π(ω). This is speciﬁed with a single parameter ν0 inverted gamma density. We use two
diﬀerent ν0 values: 25 and 5. ν0 = 5 should be interpreted as an indication of larger heteroskedasticity
compared to ν0 = 25.5 So far as the regression coeﬃcients β and the error term variance scaleparameter σ 2 are concerned, MCMC procedures are the same as homoskedastic NLR model stated in
the above paragraph.

4.2.1

Comparison of the Models: DIC and Posterior Evaluation

Since Spiegelhalter et al. (2002) ﬁrst proposed DIC it has become one of the most frequently usedmodel comparison criterion for Bayesians.6 Let a single parameter of interest be θ, then the devianceis deﬁned to be D(θ) = −2 log(likelihood). The posterior mean deviance, D(θ), is deﬁned to be D(θ)
evaluated under posterior measure, thus D(θ) = Eθ [D(θ)|y], where y denotes data. The penalty termfor DIC, pD , is deﬁned to be pD = D(θ) − D(θ̄), where θ̄ denotes posterior mean of θ. Then DIC isgiven by
DIC = D(θ) + pD .

(26)

Note that pD is often interpreted to be the eﬀective number of parameters, hence a measure of model
complexity.
The following shows D(θ), pD and DIC of the three models.
Table 1. DIC Comparedhomoskedastic NLR

Geweke (ν0 = 25)

Geweke (ν0 = 5)

Our model

D(θ)-52.975-79.049-81.598-114.956

pD

1.924

3.070

2.438

2.907

DIC

-51.051-75.980-79.160-112.050

Since the null model in (24) essentially contains two parameters, σ 2 and ϕ, pD should approximately
two. We observe that homoskedastic NLR has pD ≈ 2 while Geweke (1993)’s model and ours havepD ≈ 3. As regards to the goodness of ﬁt measure, i.e., D(θ), our model has the best measure while
5
6

See Koop (2003, p.129) for setting ν0 .
See Spiegelhalter et al. (2014) for an excellent survey, pros and cons of DIC.

–12–

the other two have lesser D(θ). In summary, our model has the best DIC with reasonable pD values
in the models compared.
We next examine three model’s parameter θ estimation result. Numerical results, except for ω
are given in Table 2. Homoskedastic NLR has the smallest ϕ and the largest σ 2 of all the models
considered. This may indicate a tradeoﬀ between ϕ and σ 2 due to misspeciﬁcation of homoskedasticity,
since other models’ ϕ and σ 2 are more or less the same. Let estimated parameter value over its standard
deviation be “t-value.” We observe that our model’s two parameters presented in Table 2, have thelargest t-value. This fact may imply GLS-like treatment of our model is most appropriate for thedataset.
Table 2. Posterior Mean Comparedhomoskedastic NLR

Geweke (ν0 = 25)

Geweke (ν0 = 5)

Our model

0.111

0.138

0.141

0.118

ϕ

(0.026)(0.025)(0.024)(0.013)

σ2

0.021

0.010

0.008

0.014

(0.004)(0.004)(0.003)(0.003)

−−

4.967

3.349

ν

(5.174)(1.735)−−“(·)” indicates posterior standard deviation.

4.3

ω: Posterior Mean and Posterior Density

We may now present 50 posterior mean’s of our remaining parameter, ω, i.e., E (ω |data) ∼ 50 × 1.
θ
Since we do not know the true heteroskedasticity of our data, however, presenting 50 posterior means
would not contribute to our understanding of unknown skedastic structure.
“Size eﬀect” hypothesis that indicates an inverse relationship between volatility to size of a company, may be an appropriate hypothesis to be dealt with. In the U.S. and world wide, starting witha seminal paper by Banz (1981) many observed “size eﬀect” that is the smaller the company is thehigher its return7 . This phenomenon could be naively associated to the mean-variance eﬃciency toyield 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. Berk (1997) among others examined the so-called size
eﬀect and proposed that a size of a company should not be measured by the market value of its equity
7

For the size eﬀect in Japan, see e.g., Chan and Chen (1991), among others.

–13–

Figure 1: Volatility versus Size: σ 2 ωi ’s and MVE(MVE) but some other variables such as sales. 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 (in one billionyen)8 . We then drew a graph with the MVE on the horizontal axis, and posterior mean of σ 2 ωi ’sfrom our model on the vertical axis. This is shown in Figure 1. The downward sloping solid line in
Figure 1, is the OLS estimated line. This ﬁgure clearly shows the larger the MVE, the smaller thevolatility as measured by σ 2 ωi . In summary, we have conﬁrmed that E (ω |data) from our model,
θgives reasonable values.
We need to see each simulated posterior density of ωi has a shape that is reasonable as a densityof variance, e.g., gamma density. To this end we selected two stocks (i) that has large MVE and small
σ 2 ωi , and (ii) that has small MVE and large σ 2 ωi , to see what the marginal posterior pdf’s of σ 2 ωiof these companies look like. Speciﬁcally, we chose Taisho Pharmaceutical Co., Ltd. for (i), and Site
Support Institute Co., Ltd. for (ii). They are given in Figures 2 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

8

Berk (1997) among others, examined the so-called size eﬀect and proposed that a size of a company should not bemeasured by MVE but some other variables such as sales.

–14–Figure 2: Marginal Posterior pdf of σ 2 ωi for Taisho Phamaceutical Co. and Site Support Institutereasonable shapes. We conclude that our Bayesian estimation of volatility supports the view that thesmaller the size of the stock, the larger is the volatility.

5

Concluding Remarks

In this paper, we proposed a Bayesian method to estimate regression error heteroskedasticity structurethat is unknown. “Unknown” in the sense that no structure is assumed. Geweke (1993)’s model issuch that when Gamma priors with a particular set of hyper parameter values are assumed on the
precision parameters of heteroskedastic regression error term, then this leads to a homoskedastic
Student-t regression error term. We pointed out that assuming such priors, are in eﬀect, imposing an
unwanted structure in the heteroskedasticity. We have, thus, proposed to use a Dirichlet prior withequal hyper parameter values. This should represent we “know nothing” status about the structureof heteroskedasticity.
We, on the other hand, believe that the Eicker-White HCCM should provide valuable information
about the heteroskedasticity, although derived from a sampling theory point of view. In empiricalanalysis, regression equation is bound to be misspeciﬁed. HCCM, in essence, draws heteroskedasticity
information connecting with regressors. Our idea is to use this HCCM information in the proposal
distribution in the Independence sampler. We showed that this approach is reasonably successful.

–15–

Finally we compared homoskedastic NLR, Geweke (1993)’s model, and our model in terms of DIC,
posterior mean signiﬁcance (we used “t-value” to this eﬀect), and posterior pdf’s. All the exercises
indicate that our model is deﬁnitely capable of drawing unknown heteroskedasticity structure.

References[1] Amemiya, T., 1985, Advanced Econometrics, Harvard University Press.
[2] Banz, R.W., 1981, The relationship between return and market value of common stocks. Journalof Financial Economics, 9, 3–18.
[3] Berk, J.B., 1997, Does size really matter? Financial Analysts Journal, 53, 12–18.
[4] Chan, K.C. and Nai-Fu Chen, 1991, Structural and return characteristics of small and large ﬁrms.
Journal of Finance, 46, 1467–84.
[5] Cochrane, J.H., 2001, Asset Pricing, Princeton University Press.
[6] Eicker, F., 1963, Asymptotic normality and consistency of the least squares estimators for familiesof linear regressions. Annals of the Mathematical Statistics, 34, 447–456.
[7] Greene, W.H., 2012, Econometric Analysis, 7th ed., Pearson.
[8] Greenberg, E., 2013, Introduction to Bayesian Econometrics, 2nd ed. Cambridge University Press.
[9] Geweke, J., 1993, “Bayesian treatment of the independent Student-t linear model,” Journal ofApplied Econometrics, 8, Supplement (December), S19–S40.
[10] Geweke, J., 2005, Contemporary Bayesian Econometrics and Statistics, Wiley.
[11] Godfrey, L.G., 2006, Tests for regression models with heteroskedasticity of unknown form. Computational Statistics & Data Analysis, 50, 2715–2733.
[12] Koop, G., 2003, Bayesian Econometrics, Wiley.
[13] Lancaster, T., 2004, An Introduction to Modern Bayesian Econometrics, Blackwell.
[14] Long, J.S. and L.H. Ervin, 2000, Using heteroscedasticity consistent standard errors in the linear
regression model. American Statistician, 54, 217–224.
[15] Robinson, P.M., 1987, Asymptotically eﬃcient estimation in the presence of heteroskedasticity of
unknown form. Econometrica, 55, 875–891.

–16–[16] Spiegelhalter, D.J., Best, N.G., Carlin, B.P. and van der Linde, A., 2002, Bayesian measures ofmodel complexity and ﬁt (with discussion), Journal of the Royal Statistical Society Series B., 64,
583–639.
[17] Spiegelhalter, D.J., Best, N.G., Carlin, B.P. and van der Linde, A., 2014, Deviance information
measure: after 12 years (with discussion), Journal of the Royal Statistical Society Series B., 76,part 3, 485–493.
[18] White, H., 1980, A heteroskedasticity-consistent covariance matrix estimator and a direct test forheteroskedasticity. Econometrica, 48, 817–838.

–17–