Gibbs estimation of microstructure models: Teaching notes

Joel Hasbrouck

Text  @ Style["Evaluated "<>DateString["DateTime"],"Subsubtitle"]

. Introduction

This note discusses Gibbs estimation of the Roll model and various modifications. The goal is a more discursive and heuristic treatment of material covered in Hasbrouck (2009). Other applications of Gibbs samplers in market microstructure include Hasbrouck (1999) and Ball and Chordia (2001).

The techniques discussed here follow an approach that relies on simulation to characterize model parameters. Applied to microstructure models, there are three key elements:

Bayesian analysis

Simulation

Characterization of microstructure data generating processes by their conditional probabilities.

Specifically:

Bayesian analysis

The models are estimated in a Bayesian framework. The differences between Bayesian and classical analysis are continually debated and discussed. The framework here is definitely Bayesian, but it should be noted that even if one doesn't buy the full Bayesian philosophy, the techniques discussed here can be motivated on grounds of estimation simplicity and computational efficiency.
This is an unusual claim. Bayesian analyses are usually more complex (both conceptually and computationally) than their classical counterparts. This is sometimes cited by Bayesian adherents as the prime barrier to their widespread adoption. Most microstructure models, though, are dynamic (over time) and they have latent (hidden, unobservable) quantities. The classic Roll model is dynamic, and the trade direction indicator ("buy or sell") variables are not observed.

Dynamic latent variable models can be formulated in state-space form and estimated via maximum likelihood. The latent variables are often non-Gaussian (e.g., again, the trade direction indicator variables), and if one wants to go beyond the techniques of multivariate linear models (like VARs), estimation involves nonlinear filtering. The Gibbs estimates are usually quicker and simpler.

There are presently a number of Bayesian statistics textbooks available. In my opinion the most useful for financial econometricians, are those that discuss econometrics from a Bayesian perspective. Lancaster (2004) and Geweke (2005) are both excellent. Lancaster's treatment is particularly accessible; Geweke presents more results. Nelson and Kim (2000) is a good introduction to the techniques in the context of a specific problem (regime switching models). In financial econometrics, the heaviest use of Bayesian simulation has been in modeling stochastic volatility. Shephard (2005) is a good survery of this area. Tanner (1996) and Carlin and Louis (2004) consider a broader range of Bayesian statistical tools and applications.

Simulation.

The output of a classical procedure (e.g., OLS) is usually a statement about the distribution of a parameter. E.g., "θ is asympototically normally distributed with mean (OverBar[θ]) and variance (Subscript[σ, θ]^2)," where the mean and variance quantities are computed directly. But we could also characterize a distribution by a sample of draws from that distribution. This is what most of modern Bayesian analysis does. The output of the estimation procedures discussed here is a stream of random draws of the parameters of interest (conditional on the model and the data). From this stream we can construct an estimate of the full distribution (via kernel smoothing) or simply a summary measure (like the mean or median).

Among other things, simulation facilitates characterization of distributions for functions of random variables.  For example, suppose that x~N(μ,(σ^2)) and we'd like to characterize the distribution of y=f(x) where f is sufficiently complicated that we can't get closed-form results.  We simply generate random values (Subscript[x, i]) and empirically examine the distribution of (Subscript[y, i])=f((Subscript[x, i])).

The link between simulation and Bayesian analysis is strong for the following reason. The distributions that arise in Bayesian analysis often describe many random variables (i.e., they are of high dimension). It also often happens that they have no closed form representation. Instead, they are characterized by simulation. The Gibbs procedure belongs to a class of random number generators called Markov Chain Monte Carlo (MCMC) techniques. They work by setting up rules for moving from one realization (draw) of the random variables to a subsequent realization. These draws are viewed as "states" in a Markov chain, and the rules define the transition probabilities. The limiting distribution of the states is identical to the distribution of the variables of interest, and is approximated by repeated application of the transition function.

Conditional probabilities

To set up a Gibbs estimate, we need to compute conditional densities for all of the unknowns (parameters and latent data). The conditional distributions for the parameters are usually identical to those found in many other applications (e.g., the normal Bayesian linear regression model). This note merely summarizes these distributions, refering the reader elsewhere for a fuller treatment. The conditional distributions for the latent data, though, are specific to the problem at hand. Although not particularly complicated, they are non-standard, and this note covers them in greater detail.

Programs

This note is written in Mathematica. Some Mathematica code and results are embedded. Most of the results, though, are computed using SAS code that is available on my web site.  These programs are available in my ftp directory at http://pages.stern.nyu.edu/~jhasbrou/Teaching/2010%20 PhD %20 Microstructure/PhDMicroSp2010/.

The programs make heavy use of SAS/IML ("Interactive Matrix Language"). This is not the language I've used for most of my papers, but it is widely available. Anyone who has a copy of SAS should be able to run the programs. These programs are not "industrial strength". I've played around with them in generating the results for this note, but they haven't been tested against all the things that might come up in, say, the CRSP daily file. I haven't done any performance benchmarks, but I suspect that they run slower than comparable code in OX or Matlab.

The main programs used here are:

RollGibbs2Trade case.sas
RollGibbsAnalyzeq.sas
RollGibbs01.sas
RollGibbsBeta01.sas

These programs call two macros: RollGibbs.sas and RollGibbsBeta.sas. These macros, in turn, make use of IML subroutines contained in a library called "imlstor". To set up this library, run the program RollGibbsLibrary01.sas (which contains the code for the subroutines).

. Overview

This note illustrates the estimation approach for the Roll (1984) model of transaction prices. In this model, the "efficient price" ((Subscript[m, t])) is assumed to follow a Gaussian random walk:

(Subscript[m, t])=(Subscript[m, t-1])+(Subscript[u, t]) where (Subscript[u, t])~N(0,(Subscript[σ, u]^2))

The transaction price ((Subscript[p, t])) is the efficient price, plus or minus a cost component that depends on whether the customer is buying or selling:

(Subscript[p, t])=(Subscript[m, t])+c (Subscript[q, t])

where c is the cost parameter and (Subscript[q, t])=±1. (If the customer is buying, (Subscript[q, t])=+1; if selling, (Subscript[q, t])=-1). The trade prices are observed. The (Subscript[q, t]) and the (Subscript[m, t]) are not. Taking first differences:

Δ (Subscript[p, t])=c (Subscript[Δq, t])+(Subscript[u, t])

This specification is important because if the (Subscript[Δq, t]) were known, this would be a simple regression.

Bayesian estimation of normal linear regressions is well understood. The discussion (in the next section) starts with a review of these procedures.

There are two parameters in this "regression": c (the coefficient) and (Subscript[σ, u]^2). It is fairly easy to compute (in closed form) the posterior distributions f(c|(Subscript[σ, u]^2), (Subscript[p, 1]),...,(Subscript[p, T])) and f((Subscript[σ, u]^2)|c, (Subscript[p, 1]),...,(Subscript[p, T])).  It is not possible to compute in closed form the joint posterior f(c,(Subscript[σ, u]^2)| (Subscript[p, 1]),...,(Subscript[p, T])). This motivates the next section, which summarizes the Gibbs sampler.

The Gibbs procedure is illustrated by applying it to a special case of the Roll model, one in which c and (Subscript[σ, u]^2) are known, but the (Subscript[q, t]) are not. The note then turns to a full estimation of the Roll model, and extensions.

. Mathematica initializations

SetDirectory[NotebookDirectory[]];

<<MVN.m

<<Notation`

The following commands define symbolizations that are convenient for labeling things.

Symbolize[(Subscript[Anything_, Rule])];Symbolize[(Subscript[Anything_, Rules])]

Off[General::"spell",General::"spell1"];

. Bayesian analysis of the normal linear regression model

The basic Bayesian approach

Bayesian analysis consists of using a model and data to update prior beliefs.The revised beliefs are usually called posterior beliefs, or simply "the posterior". Let y denote the observed data, and let the model be specified up to a parameter θ (possibly a vector). Bayes theorem says :

f(θ|y)=((f(θ,y))/(f(y)))=((f(y|θ)f(θ))/(f(y)))∝f(y|θ)f(θ)

f(θ) is an assumed prior distribution for the parameter.

f(y|θ) is the likelihood function for the observations, given a particular value of θ.

The use of ∝ ("is proportional to") reflects the fact that it is usually not necessary to compute f(y), at least not by computing the marginal f(y)=∫f(θ,y)d θ. Instead, f(y) can treated as a normalization constant, set so that the posterior integrates to unity.

Often a distribution of interest, say f(x),  can be written as f(x)=k g(x), where g(x) is a parsimonious function of x and k is some scaling factor. k might in fact be very complicated, possibly depending on other random variables and implicitly incorporating other distribution functions, but for purposes of characterizing the distribution of x, it is constant. In this case, g(x) is said to be the kernel of f(x).

Bayesian estimation of the normal linear regression model

The normal regression model is:

(Underscript[y, N×1])=(Underscript[X, N×K]) (Underscript[β, K×1]) + (Underscript[u, N×1]) where u~N(0,(Subscript[Ω, u]))

X is a matrix of covariates (explanatory variables) possibly including a constant; β is the coefficient vector.

Estimation of coefficients (given the error variance)

Assume for the moment that (Subscript[σ, u]^2) is known. It is particularly convenient to assume a multivariate normal prior distribution for the coefficients:

β~N((Subscript[μ, β]^Prior),(Subscript[Ω, β]^Prior))

The posterior distribution, f(β|y)is

N((Subscript[μ, β]^Post),(Subscript[Ω, β]^Post))

where (Subscript[μ, β]^Post)=D d; 	(Subscript[Ω, β]^Post)=(D^(-1));  	(D^(-1))=X'(Subscript[Ω, u]^(-1))X+((((Subscript[Ω, β]^Prior)))^(-1));  	d=X'(Subscript[Ω, u]^(-1))y+((((Subscript[Ω, β]^Prior)))^(-1))(Subscript[μ, β]^Prior)

As (Subscript[Ω, β]^Prior) increases in magnitude, the posterior mean and variance converge toward the usual OLS values.
In this case, both the prior and posterior have the same form (multivariate normal). Such a prior is said to be conjugate.

Simulating the coefficients

We'll often have to make a random draw from the coefficient distribution. To make a random draw from 
(Underscript[x, n×1])~MVN(μ,Ω):

Compute the Cholesky factorization F: Ω=F'F, where F is an upper triangular matrix.

Draw (Underscript[z, n×1]) where the (Subscript[z, i]) are i.i.d. N(0,1).

Set the random draw as x=μ+F'z

Restrictions on the prior

The economic model sometimes imposes bounds on the coefficients. For example, in the Roll model, we'll often want to require c>0. Suppose that the coefficient prior is

β~N((Subscript[μ, β]^Prior),(Subscript[Ω, β]^Prior)), (UnderBar[β])<β<(OverBar[β])

Note that when we write this, (Subscript[μ, β]^Prior),(Subscript[Ω, β]^Prior) denote the formal parameters of the normal density. But since the distribution is truncated, they no longer denote the mean and covariance of the density.

With this prior, the posterior is simply N((Subscript[μ, β]^Post),(Subscript[Ω, β]^Post)), with (Subscript[μ, β]^Post) and (Subscript[Ω, β]^Post)computed as described above, restricted to the space (UnderBar[β])<β<(OverBar[β]).

Simulation from restricted normals.

First suppose that we want to make a random draw z from a standard normal density, restricted to the interval (UnderBar[z])<z<(OverBar[z]). The procedure is:

Compute (UnderBar[π])=Φ((UnderBar[z])) and (OverBar[π])=Φ((OverBar[z])), where Φ is the c.d.f. of the standard normal.

Draw u from the uniform distribution over ((UnderBar[π]),(OverBar[π]))

Set z=(Φ^(-1))(u).

Now suppose that we want to make a bivariate random draw from x=(({{(Subscript[x, 1])}, {(Subscript[x, 2])}}))~N(μ,Ω), (UnderBar[x])<x<(OverBar[x]).

Compute the Cholesky factorization F: Ω=F'F, where F is an upper triangular matrix.

Set (UnderBar[z])=(((UnderBar[(Subscript[x, 1])])-(Subscript[μ, 1]))/(Subscript[F, 11])) and (OverBar[z])=(((OverBar[(Subscript[x, 1])])-(Subscript[μ, 1]))/(Subscript[F, 11]))

Draw (Subscript[z, 1]) from the standard normal density, restricted to ((UnderBar[z]),(OverBar[z])). Then (Subscript[x, 1])=(Subscript[μ, 1])+(Subscript[F, 11])(Subscript[z, 1]) will have the properties required of (Subscript[x, 1]).

Set η=(Subscript[F, 11])(Subscript[z, 1]).

Set (UnderBar[z])=(((UnderBar[(Subscript[x, 2])])-(Subscript[μ, 2])-η)/(Subscript[F, 22])) and (OverBar[z])=(((OverBar[(Subscript[x, 2])])-(Subscript[μ, 2])-η)/(Subscript[F, 22]))

Draw (Subscript[z, 2]) from the standard normal density, restricted to ((UnderBar[z]),(OverBar[z])). Then (Subscript[x, 2])=(Subscript[μ, 2])+(Subscript[F, 22])(Subscript[z, 2]) will*********** have the properties required of (Subscript[x, 2]).

The random draw as x=μ+F'z will have the required joint properties

This method may be generalized to higher dimensions. (See Hajivassiliou, V., D. McFadden and P. Ruud (1996). Simulation of multivariate normal rectangle probabilities and their derivatives - Theoretical and computational results. Journal of Econometrics 72(1-2): 85-134)

Estimation of error variance (given the coefficients)

Assuming that β is known, it is convenient to specify a inverted gamma prior for (Subscript[σ, u]^2). One way of writing this is:

(1/(Subscript[σ, u]^2))~Γ[(a^Prior),(b^Prior)]

Then the posterior is

(1/(Subscript[σ, u]^2))|y~Γ[(a^Post),(b^Post)]

where  (a^Post)=(a^Prior)+(N/2) and (b^Post)=(b^Prior)+(((∑_(i=1)^N)(Subscript[u, i]^2))/2)

The (Subscript[u, i]) are the regression residuals u=y-X β.

Further notes

The density of a gamma variate x with parameters a and b is:

PDF[GammaDistribution[a,λ],x]/.λ->1/b

Mean[GammaDistribution[a,λ]]/.λ->1/b

Note: In the statistics literature, the Gamma distribution with parameters aand b is usually expressed as immediately above. Mathematica parameterizes the distribution with the second parameter expressed as an inverse.

PDF[GammaDistribution[a,λ],x]/.λ->1/b/.x->1/z

. The Gibbs recipe

The Gibbs procedure is a vehicle for simulating from a complicated joint distribution f((Subscript[x, 1]),...,(Subscript[x, n])), possibly one that possesses no closed form representation.

The draws are constructed by iterating over the full conditional distributions:

f((Subscript[x, 1])|(Subscript[x, 2]),...,(Subscript[x, n]))
f((Subscript[x, 2])|(Subscript[x, 1]),(Subscript[x, 3]),...,(Subscript[x, n]))
...
f((Subscript[x, n])|(Subscript[x, 1]),...,(Subscript[x, n-1]))

Let x=(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {⋮}, {(Subscript[x, n])}}))
Each iteration of the Gibbs sampler is called a sweep. 
Let (x^([i])) denote the value of x at the conclusion of the ith sweep.
The procedure is:

Initialization. Set (x^([0])) to any feasible value.

Sweep i:
Given (x^([i-1])):

Draw (Subscript[x, 1]^([i])) from f((Subscript[x, 1])|(Subscript[x, 2]^([i-1])),...,(Subscript[x, n]^([i-1])))

Draw (Subscript[x, 2]^([i])) from f((Subscript[x, 2])|(Subscript[x, 1]^([i])),(Subscript[x, 3]^([i-1])),...,(Subscript[x, n]^([i-1])))
...

Draw (Subscript[x, n]^([i])) from f((Subscript[x, n])|(Subscript[x, 1]^([i])),...,(Subscript[x, n-1]^([i])))

Repeat

In the limit, as i->∞,  (x^([i]))is distributed as f(x).

Notes

The (x^([i])) are not independently distributed: (x^([i])) takes (x^([i-1])) as its starting point. If the degree of dependence is high, a large number of sweeps may be needed to ensure proper mixing.

The dependence may affect the calculation of some summary statistics. Think of the analogy to standard time series analysis. If (Subscript[z, 1]),...,(Subscript[z, T]) are a sample of stationary time series data, (Subscript[Σz, i])/T is a consistent estimate of E[(Subscript[z, t])]. The standard error of this estimate, however, must be corrected for dependence.

Convergence may be an issue. It is useful to graph the full sequence of draws.

In analyzing the sequence of draws, it is common to throw out a few initial draws, so as to reduce the dependence on starting values. These discarded draws are sometimes called burn in draws.

The Gibbs sampler also works when multiple variables are drawn at once. We might, for example, draw (Subscript[x, 1]^([i])) and (Subscript[x, 2]^([i])) from f((Subscript[x, 1]),(Subscript[x, 2])|(Subscript[x, 3]^([i-1])),...,(Subscript[x, n]^([i-1]))). This block sampling is often more computationally efficient.

Application to the normal regression model

From earlier results, we have f(β|y,(Subscript[σ, u]^2)) and f((Subscript[σ, u]^2)|β,y). To obtain the full posterior f(β,(Subscript[σ, u]^2)|y) via the Gibbs procedure:

Initialize (Subscript[σ, u]^(2[0])) to any feasible value. The ith sweep of the sampler is:

Draw (β^([i])) from f(β|y,(Subscript[σ, u]^(2[i-1]))). (This will be a draw from a multivariate normal posterior.)

Draw (Subscript[σ, u]^(2[i])) from f[(Subscript[σ, u]^2)|y,(β^([i]))). (That is, draw 1/(Subscript[σ, u]^(2[i])) from the gamma posterior.)

Proceed, iteratively drawing (Subscript[σ, u]^(2[i])) and (β^([i])).

Notes

The f(β,(Subscript[σ, u]^2)|y) is an exact small-sample distribution.

We now return to ...

. The Roll model

Recall that we're working with the price change specification:

Δ (Subscript[p, t])=c Δ (Subscript[q, t])+(Subscript[u, t])

The sample is (Subscript[p, 1]),(Subscript[p, 2]),...,(Subscript[p, T]), and there are T-1 price changes.

The unknowns are the parameters c and (Subscript[σ, u]^2), and the latent data: (Subscript[q, 1]),...,(Subscript[q, T]).
In the Bayesian perspective, parameters and latent data are treated identically, and "estimated" in a similar fashion.
We don't need to construct priors for the (Subscript[q, t]). We can use the ones we assumed in the data generating process: (Subscript[q, t])=±1 with equal probabilities.
The prior on c is c~N((μ^Prior),(Ω^Prior)) restricted to c>0. (I often take (μ^Prior)=0 and (Ω^Prior)=1. Remember that these are parameters of the truncated distribution, not the true mean and variance.)
The prior on (Subscript[σ, u]^2) is inverted gamma with parameters a and b. (I often take a=b=(10^(-6)).)
The Gibbs sampler will look like this:

Initialize (c^([0])), (Subscript[σ, u]^(2[0])) and (Subscript[q, 1]^([0])),..., (Subscript[q, T]^([0])) to any feasible values.

(I usually take (Subscript[q, 1])=1; (Subscript[q, t])=Sign((Subscript[Δp, t])) if (Subscript[Δp, t])≠0, (Subscript[q, t])=(Subscript[q, t-1]) if (Subscript[Δp, t])=0. For US equities, c=0.01 ("1%") is a good ballpark figure, if we're working in logs, and (Subscript[σ, u]^2)=(0.01^2)).

For the ith sweep of the sampler:

Estimate the price change specification as a regression, assuming that (Subscript[q, t])=(Subscript[q, t]^([i-1])) and that (Subscript[σ, u]^2)=(Subscript[σ, u]^(2[i-1])). 
Construct the posterior for c, and draw (c^([i])) from this posterior.

Using (c^([i])), compute the residuals from the regression. Construct the posterior for (Subscript[σ, u]^2) and draw (Subscript[σ, u]^(2[i])) from this posterior.

Draw (Subscript[q, 1]^([i])) from f((Subscript[q, 1])|(c^([i])),(Subscript[σ, u]^(2[i])),(Subscript[q, 2]^([i-1])),(Subscript[q, 3]^([i-1])),..., (Subscript[q, T]^([i-1])))
Draw (Subscript[q, 2]^([i])) from f((Subscript[q, 2])|(c^([i])),(Subscript[σ, u]^(2[i])),(Subscript[q, 1]^([i])),(Subscript[q, 3]^([i-1])),..., (Subscript[q, T]^([i-1])))
...
Draw (Subscript[q, T]^([i])) from f((Subscript[q, T])|(c^([i])),(Subscript[σ, u]^(2[i])),(Subscript[q, 1]^([i])),(Subscript[q, 2]^([i])),..., (Subscript[q, T-1]^([i])))

The first two steps of the sweep are, as discussed, standard Bayesian procedures. We now turn to the third step: simulating the trade direction indicators

. Estimating the (Subscript[q, t]) using a Gibbs sampler.

In this section, we'll be taking c and (Subscript[σ, u]^2) as known. We'll first look at the simple case where T=2. We can get closed-form results here, so we don't really need a Gibbs procedure, but it's a good starting point.

The distribution of (Subscript[q, 1]) and (Subscript[q, 2]) when T=2

Symbolize[(Subscript[p, _])]; Symbolize[(Subscript[q, 1])];Symbolize[(Subscript[q, 2])]; Symbolize[(Subscript[σ, _])]; Symbolize[(Subscript[u, 2])];

$TextStyle={FontFamily->"Times",FontSize->12,LineSpacing->{1,-1}};

Graph

s=0.05; hc=0;hl=0.02;

Suppose:

(Subscript[p, 1])=1; (Subscript[p, 2])=1.8; c=.5;

ListPlot[{{1,(Subscript[p, 1])+c},{1,(Subscript[p, 1])-c},{2,(Subscript[p, 2])+c},{2,(Subscript[p, 2])-c}},PlotRange->{{0.5,2.5},{0,2.5}},
AxesOrigin->{0.5,0},PlotStyle->PointSize[0.01],
Frame->True,FrameLabel->{Time,Price,None,None},FrameTicks->{{0,1,2},Automatic,None,None},
Background->White,
Epilog->{
Text["(Subscript[p, 1])",{1-s,(Subscript[p, 1])},{1,0}],
Text["(Subscript[m, 1])((Subscript[q, 1])=-1)",{1-s,(Subscript[p, 1])+c},{1,0}],
Text["(Subscript[m, 1])((Subscript[q, 1])=+1)",{1-s,(Subscript[p, 1])-c},{1,0}],
Text["(Subscript[p, 2])",{2+s,(Subscript[p, 2])},{-1,0}],
Text["(Subscript[m, 2])((Subscript[q, 2])=-1)",{2+s,(Subscript[p, 2])+c},{-1,0}],
Text["(Subscript[m, 2])((Subscript[q, 2])=+1)",{2+s,(Subscript[p, 2])-c},{-1,0}],Point[{1,(Subscript[p, 1])}],Point[{2,(Subscript[p, 2])}]}]

ListPlot[{{1,(Subscript[p, 1])+c},{1,(Subscript[p, 1])-c},{2,(Subscript[p, 2])+c},{2,(Subscript[p, 2])-c}},PlotRange->{{0.5,2.5},{0,2.5}},
AxesOrigin->{0.5,0},PlotStyle->PointSize[0.01],
Frame->True,FrameLabel->{Time,Price,None,None},FrameTicks->{{0,1,2},Automatic,None,None},
Background->White,
Epilog->{
Arrowheads[Medium],
Arrow[{{1,(Subscript[p, 1])+c},{2,(Subscript[p, 2])+c}},.05],
Arrow[{{1,(Subscript[p, 1])+c},{2,(Subscript[p, 2])-c}},.05],
Arrow[{{1,(Subscript[p, 1])-c},{2,(Subscript[p, 2])+c}},.05],
Arrow[{{1,(Subscript[p, 1])-c},{2,(Subscript[p, 2])-c}},.05],
Text["A. (Subscript[u, 2])((Subscript[q, 1])=-1, (Subscript[q, 2])=-1)",{1.5,2},{1,0}],
Text["B. (Subscript[u, 2])((Subscript[q, 1])=-1, (Subscript[q, 2])=+1)",{1.7,1.5},{-1,0}],
Text["C. (Subscript[u, 2])((Subscript[q, 1])=+1, (Subscript[q, 2])=-1)",{1.9,2},{-1,0}],
Text["D. (Subscript[u, 2])((Subscript[q, 1])=+1, (Subscript[q, 2])=+1)",{1.5,.7},{-1,0}]}]

Given (Subscript[p, 1]) and c, a choice of (Subscript[q, 1]) determines m, and similarly for (Subscript[p, 2]). Therefore, setting (Subscript[q, 1]) and (Subscript[q, 2]) fixes (Subscript[u, 2])=(Subscript[m, 2])-(Subscript[m, 1]).

Since (Subscript[q, 1]),(Subscript[q, 2])∈{-1,+1}, there are four possible paths (values) for (Subscript[u, 2]). These are labeled A, B, C and D in the figure. Since low values of |(Subscript[u, 2])| are more likely than high values, the relative lengths of these paths indicate the relative likelihood of the ((Subscript[q, 1]),(Subscript[q, 2])) realizations that determined them:

Intuitively, path B is the shortest, so it is most likely that (Subscript[q, 1])=-1, (Subscript[q, 2])=+1.

Path C is the longest, so it is least likely that (Subscript[q, 1])=+1, (Subscript[q, 2])=-1.

Paths A and D are of equal length, corresponding to the realizations (Subscript[q, 1])=(Subscript[q, 2])=-1 and (Subscript[q, 1])=(Subscript[q, 2])=+1.

We now turn to a more exact treatment.

Clear[s,(Subscript[p, 1]),(Subscript[p, 2]),c];

The joint density of (Subscript[q, 1]) and (Subscript[q, 2])

The density function for u is f((Subscript[u, 2]))=

PDF[NormalDistribution[0,(Subscript[σ, u])],u]

By rearranging the price change specification, (Subscript[u, 2])=

(Subscript[u, 2])[(Subscript[q, 1]):_,(Subscript[q, 2]):_]:=(Subscript[Δp, 2])-c((Subscript[q, 2])-(Subscript[q, 1]))

(Subscript[u, 2])[(Subscript[q, 1]),(Subscript[q, 2])]

The probability Pr[(Subscript[q, 1]),(Subscript[q, 2])]∝f((Subscript[u, 2])[(Subscript[q, 1]),(Subscript[q, 2])])

(Subscript[Pr, Rule])=Pr[(Subscript[q, 1]):_,(Subscript[q, 2]):_]:>PDF[NormalDistribution[0,(Subscript[σ, u])],(Subscript[u, 2])[(Subscript[q, 1]),(Subscript[q, 2])]]//Simplify;

Pr[(Subscript[q, 1]),(Subscript[q, 2])]/.(Subscript[Pr, Rule])

The possible outcomes are:

Outcomes={{-1,-1},{-1,+1},{+1,-1},{+1,+1}}; 
TableForm[Outcomes,TableHeadings->{None,{(Subscript[q, 1]),(Subscript[q, 2])}}]

We'll normalize by the sum of the probabilities:

PrSum=Plus @@ Apply[Pr,Outcomes,{1}]

(Subscript[PrN, Rule])=PrN[(Subscript[q, 1]):_,(Subscript[q, 2]):_]:>Pr[(Subscript[q, 1]),(Subscript[q, 2])]/PrSum;

The normalized probability is:

PrN[(Subscript[q, 1]),(Subscript[q, 2])]/.(Subscript[PrN, Rule])

For demonstration purposes, here are some values:

nValues={(Subscript[Δp, 2])->.8,(Subscript[σ, u])->1,c->0.5};

With thes values, the normalized probabilities are:

TableForm[Transpose[Append[Transpose[Outcomes],Apply[PrN,Outcomes,{1}]/.(Subscript[PrN, Rule])/.(Subscript[Pr, Rule])/.nValues]],
TableHeadings->{None,{(Subscript[q, 1]),(Subscript[q, 2]),"Probability"}}]

Gibbs sampler

We can draw the trade direction indicator variables jointly in this case. There's no need to use a Gibbs sampler. But for illustration purposes, let's build it.

The required conditional probabilities are Pr[(Subscript[q, 1])|(Subscript[q, 2])] and Pr[(Subscript[q, 2])|(Subscript[q, 1])].

These may be computed directly from joint distribution given above, but it is usually computationally easier to work with the odds ratio.

For example, (Subscript[u, 2])[(Subscript[q, 1]),(Subscript[q, 2])]=

(Subscript[u, 2])[(Subscript[q, 1]),(Subscript[q, 2])]

Its density is:

f=PDF[NormalDistribution[0,(Subscript[σ, u])],(Subscript[u, 2])[(Subscript[q, 1]),(Subscript[q, 2])]]

So given (Subscript[q, 2]), the odds in favor of a buy at time 1 are Odds(Buy)=((Pr[(Subscript[q, 1])=+1|...])/(Pr[(Subscript[q, 1])=-1|...]))=

OddsBuy1=((f/.(Subscript[q, 1])->+1)/(f/.(Subscript[q, 1])->-1))//Simplify

Then Pr[Buy]=((Odds(Buy))/(1+Odds(Buy))). We compute this probability and make a draw for (Subscript[q, 1]).

For the particular numeric values we worked with above ((Subscript[Δp, 2])=0.8, c=0.5, (Subscript[σ, u])=1), these odds, for (Subscript[q, 2])=+1 and (Subscript[q, 2])=-1 are:

OddsBuy1/.(Subscript[Δp, 2])->0.8/.c->0.5/.(Subscript[σ, u])->1/.(Subscript[q, 2])->{+1,-1}

So, for example, if (Subscript[q, 2])=+1, Pr[(Subscript[q, 1])=+1|...]=

%[[1]]/(1+%[[1]])

Similarly, given (Subscript[q, 1]), the odds in favor of a buy at time 2 are ((Pr[(Subscript[q, 2])=+1|...])/(Pr[(Subscript[q, 2])=-1|...]))=

OddsBuy2=((f/.(Subscript[q, 2])->+1)/(f/.(Subscript[q, 2])->-1))//Simplify

For the numeric values, these are (for (Subscript[q, 1])=+1 and (Subscript[q, 1])=-1):

OddsBuy2/.(Subscript[Δp, 2])->0.8/.c->0.5/.(Subscript[σ, u])->1/.(Subscript[q, 1])->{+1,-1}

The Gibbs sampler involves the following steps:

We construct a series of realizations in the following fashion:

Initialize (Subscript[q, 1]^([0])),(Subscript[q, 2]^([0]))=±1 (it doesn't matter which).

The ith sweep involves the following steps:

Draw (Subscript[q, 1]^([i])) from Pr[(Subscript[q, 1])|(Subscript[q, 2]^([i-1]))]

Draw (Subscript[q, 2]^([i])) from Pr[(Subscript[q, 2])|(Subscript[q, 1]^([i]))]

After N sweeps we'll  have a series of N simulated realizations: (q^([0])),(q^([1])),...,(q^([N])) where (q^([i]))=((Subscript[q, 1]^([i])),(Subscript[q, 2]^([i]))).
In the limit, as N->∞, the distribution of (q^([N])) is the joint distribution Pr[(Subscript[q, 1]),(Subscript[q, 2])]

The sampler was implemented in a SAS program (RollGibbs2TradeCase.sas), which was run for 10,000 sweeps. The tabulated frequencies of the simulated draws were:

Compare these to the computed probabilities above.

Why do we ever need to use Gibbs sampler when we can compute path probabilities directly? We need to compute path probabilities for the entire sample. With two price observations, there are (2^2)=4 buy/sell paths. A year contains about 250 trading days. The number of buy/sell paths is:

(2^250)

The T-trade case

A sweep of the sampler will involve:

Draw (Subscript[q, 1]^([i])) from Pr((Subscript[q, 1])|(Subscript[q, 2]^([i-1])),(Subscript[q, 3]^([i-1])),..., (Subscript[q, T]^([i-1])))

Draw (Subscript[q, 2]^([i])) from Pr((Subscript[q, 2])|(Subscript[q, 1]^([i])),(Subscript[q, 3]^([i-1])),..., (Subscript[q, T]^([i-1])))
...

Draw (Subscript[q, T]^([i])) from Pr((Subscript[q, T])|(Subscript[q, 1]^([i])),(Subscript[q, 2]^([i])),..., (Subscript[q, T-1]^([i])))

In general, Pr[(Subscript[q, t])|...]depends only on the adjacent trades - those at times t-1 and t+1. So the sampler becomes:

Draw (Subscript[q, 1]^([i])) from Pr((Subscript[q, 1])|(Subscript[q, 2]^([i-1])))

Draw (Subscript[q, 2]^([i])) from Pr((Subscript[q, 2])|(Subscript[q, 1]^([i])),(Subscript[q, 3]^([i-1])))
...

Draw (Subscript[q, T]^([i])) from Pr((Subscript[q, T])| (Subscript[q, T-1]^([i])))

The first draw, for (Subscript[q, 1]), is the same as the draw for (Subscript[q, 1]) in the T=2 case.

The last draw, for (Subscript[q, T]), is the same as the draw for (Subscript[q, 2]) in the T=2 case.

We now turn to the middle draws.

(Subscript[u, Rule]) = (Subscript[u, t_]):>(Subscript[Δp, t])+c((Subscript[q, t])-(Subscript[q, t-1]));

(Subscript[u, t+1])/.(Subscript[u, Rule])

Since the (Subscript[u, t]) are assumed to be independent, the joint density f((Subscript[u, t]),(Subscript[u, t+1]))∝

PDF[NormalDistribution[0,(Subscript[σ, u])],(Subscript[u, t])]PDF[NormalDistribution[0,(Subscript[σ, u])],(Subscript[u, t+1])]

f=%/.(Subscript[u, Rule])//Simplify

The odds ratio is:

Odds=((Pr[(Subscript[q, t])=+1|...])/(Pr[(Subscript[q, t])=-1|...]))

Odds=((f/.(Subscript[q, t])->+1)/(f/.(Subscript[q, t])->-1))//Simplify

from which we may compute

Pr[(Subscript[q, t])=+1|...]=(Odds/(1+Odds))

and make the desired draw.

Examples (SAS program RollGibbs Analyze q.sas)

I simulated twenty trades for a price process with (Subscript[σ, u])=0.01 and c=0.01, and then ran the Gibbs sampler for 2,000 sweeps to estimate the trade directions. Here is a plot of the transaction prices:

Below are the actual and average simulated trade directions. Actuals are indicated by a dot; estimated are indicated by the bars.

Note that estimates by and large agree with with actuals, at least in direction ((Subscript[q, 16]) is the sole exception).

Now, consider the same analysis, with the cost parameter changed to c=0.005, i.e., one-half the previous value. The figure below shows the prices. Notice that the bid-ask bounce is much less visually evident.

Here are the actual and estimated trade direction indicators:

The buy/sell classification accuracy here is weaker. There are more directionally incorrect inferences, and the directions that are correctly identified are weaker.  Just as we'd have a harder time picking out the buys and sells visually, the sampler has a tougher time classifying trades.

We'll now try things with (Subscript[σ, u])=0.01 and c=0.001. Here's the price path:

And here are the actual and estimated trade directions.

The intuition is as follows. Intuitively, the Gibbs sampler tries to figure out how much of an oberserved price change is transient (due to bid ask effects) and how much is permanent (the efficient price innovation). When c is large relative to (Subscript[σ, u]), bid-ask bounce generates reversals that are easy to pick out visually, and using the sampler.  When c is small, though, bid-ask effects are swamped by the innovations in the efficient price, and are not easily discerned.

We'll see that this extends to the parameter estimates as well.

Modification when some of the q's are known.

In some samples, it might happen that the trade directions are known for some subset of the (Subscript[q, t]). For these (Subscript[q, t]), we don't simulate; we simply leave them at their known values.

This might seem to violate the assumed probability structure of the model in a fundamental way. After all, if the data generating process and our priors are that (Subscript[q, t])=±1, with equal probability, how can a definite realization be accomodated? The answer is that we're conditioning on the observed data, and the only thing that matters is the prior distribution of the (Subscript[q, t]) that we don't observe.

By way of a more formal justification, we could assume that the data generating process involves two steps:

First (Subscript[q, t]) is drawn, ±1, each with probability (1/2).

Next, an indicator variable (Subscript[O, t]) is drawn. With probability π, (Subscript[O, t])=1, and the actual (Subscript[q, t]) is observed. With probability 1-π, (Subscript[O, t])=0, and the actual (Subscript[q, t]) is unobserved.

As part of the sample, we "observe" the realizations of (Subscript[O, t]). That is, we know which (Subscript[q, t]) are known for sure. If we don't care about modeling the (Subscript[O, t]) process, letting the observed (Subscript[q, t]) remain at their known values and simulating the rest corresponds to estimation conditional on the realized (Subscript[O, t]). This is a sensible way to proceed.

In doing this, we are implicitly assuming that the (Subscript[O, t]) process is independent of (Subscript[q, t]) and (Subscript[u, t]). If buys are more likely to be observed than sells, or if the realization of (Subscript[O, t]) depends on the magnitude of (Subscript[u, t]) ("Trades are more likely to be observed for large efficient price movements"), then the (Subscript[O, t]) process is informative, and we may wish to model it explicitly.

Modification when some of the q's are zero.

The U.S. daily CRSP dataset reports closing prices. But if there is no trade on day t, then the reported price is the midpoint of the closing bid and ask. This event is indicated in the data (using a negative price). Essentially, (Subscript[p, t])=(Subscript[m, t]).

Formally, we can incorporate this into the data generating process by noting that (Subscript[p, t])=(Subscript[m, t])+c (Subscript[q, t]) = (Subscript[m, t]) when (Subscript[q, t])=0.

We can handle this situation in a fashion similar to the known-(Subscript[q, t]) case discussed above. If the price reported at time t is a quote midpoint, we set (Subscript[q, t])=0, and don't resample it during the sweeps.

Formally, this can be justified by letting (Subscript[O, t]) denote an indicator variable of whether or not there was a trade. Estimation can then proceed condtional on the (Subscript[O, t]) realizations.  Here as well, we're implicitly assuming idependence. We're ruling out (or at least not modeling), for example, the possibility that trades are more likely to occur on days when there are large efficient price changes.

. Full estimation of the Roll model

Sample runs from Roll Gibbs 01.sas

In all cases, the prior on c is N((Subscript[μ, c]^Prior),(Subscript[Ω, c]^Prior)), restricted to c>0, with (Subscript[μ, c]^Prior)=0 and (Subscript[Ω, c]^Prior)=1. The prior on (Subscript[σ, u]^2) is inverted gamma with a=b=1×(10^(-6)).

(Subscript[σ, u])=.01; c=0.01, 100 observations, 20,000 sweeps (first 20% dropped)

Posteriors:

Variable        N      Mean       Dev       Min       Max
---------------------------------------------------------
SWEEP       16000     12001      4619      4001     20000
SDU         16000    0.0129    0.0016    0.0088    0.0211
C           16000    0.0091    0.0014    0.0001    0.0132
---------------------------------------------------------

Contour plot of joint posterior:

Note the downward slope of posterior. The procedure is trying to allocate volatility, either to the permanent (random-walk) component or to the transient component (effective cost). If more volatility is attributed to the sdu component, less is attributed to c.

(Subscript[σ, u])=.01; c=0.001, 100 observations, 20,000 sweeps (first 20% dropped)

When c<<(Subscript[σ, u]), the transient cost effects (reversals) are difficult to disentangle from the random-walk component. We still get a reasonably sharp posterior for (Subscript[σ, u])(=sdu), but the posterior for c is broad.

Posteriors:

Variable        N      Mean       Dev       Min       Max
---------------------------------------------------------
SWEEP       16000     12001      4619      4001     20000
SDU         16000    0.0107    0.0008    0.0082    0.0148
C           16000    0.0014    0.0010      8E-8    0.0057
---------------------------------------------------------

Contour plot of joint posterior:

Use of an unrestricted prior for c

Is it really necessary to require c>0? After all, if c really is non-negative, shouldn't the procedure pick it up?

When we run the last problem with no restrictions on c, two things happen.

For all of the (Subscript[q, t]), the sample draws average out to zero. (We can't tell whether a given trade is a buy or a sell.)

The posterior for c is bimodal, and symmetric about zero.

Here's an example. 30 observations were simulated with (Subscript[σ, u])=0.01 and c=0.01. The prior for c was not restricted to be non-negative. Here is the posterior:

What's happening is this. The Roll model with c>0 is observationally equivalent to one in which c<0, and the trade signs are reversed. Without the nonnegativity restriction on c, the posterior (and the sampler) will span both possibilities.

Using vague prior for c

In the preceding simulations, the prior for c is N(0,1), restricted to c>0. This is fairly flat over the usual region of interest. (For a US equity, even extreme values of c, estimated from trade and quote data, are rarely above 0.05.) Nevertheless, there is some curvature. Why not remove all doubt and set the prior to, say, N(0,1000000)?

The problem with this is that under some circumstances we may need to make a draw from the prior. This is not common, but in a small sample, with c<<(Subscript[σ, u]), over many draws, the following situation may arise. Suppose that on a particular sweep, the trade direction indicators are drawn to have the same sign: (Subscript[q, 1])=...=(Subscript[q, T])=+1 or (Subscript[q, 1])=...=(Subscript[q, T])=-1. In this case, all of the (Subscript[Δq, t]) are zero. But the (Subscript[Δq, t]) are the r.h.s. variables in the price change regression. If an explanatory variable in a regression has no variation, the regression is completely uninformative. In this case, the draw of c for that sweep must be made from the prior.

. Return factors

Return factors are logically introduced by adding them to the efficient price change specification, e.g.,

(Subscript[m, t])=(Subscript[m, t-1])+ (Subscript[f, t]^′)b + (Subscript[u, t])

where (Subscript[f, t]) is a K×1 vector of known factor realizations and b is the vector of loadings (factor coefficients). The factor terms then appear in the trade price specification:

Δ (Subscript[p, t])=c Δ (Subscript[q, t])+ (Subscript[f, t]^′)b + (Subscript[u, t])

For example, a market model for stock i might be specified as:

(Subscript[Δp, i t])=(Subscript[c , i])Δ (Subscript[q, it ])+ (Subscript[b, i])(Subscript[r, m]) + (Subscript[u, t])

where (Subscript[r, m]) is the return on a market index.

The b are estimated as another coefficient in the Bayesian regression (the same regression in which cis estimated). The parameters of the coefficient prior (Subscript[μ, β]^Prior) and (Subscript[Ω, β]^Prior) are expanded to include the new coefficients, and they are drawn (simulated) in the same step as the draw for c.

The addition of a return factor will generally increase the resolution between permanent (random-walk) and transient (bid-ask) components. In practice, this accomplishes two things:

We'll get a better estimate of c.

We'll also generally get better estimates of b.

The first of these is pretty straightforward: the explanatory power of the factors reduces the residual variance. The second point may require some expansion. bs are conventionally estimated using daily price changes. These daily price changes are contaminated by bid-ask bounce. For some stocks, bid-ask bounce may be large compared to the factor-induced and idiosyncratic changes in the efficient price, leading to large estimation errors in the bs.  Estimating a specification that includes a c (Subscript[Δq, t]) term effectively allows us to estimate bs on price series that are purged of the bid-ask bounce.

Market model example (RollGibbsBeta 01.sas)

The specification is the one-factor model described above The parameters are: c=0.01; (Subscript[r, m])~N(0,(Subscript[σ, m]^2)) with (Subscript[σ, m])=0.01; (Subscript[σ, u])=0.01; β=1.1.
100 observations were simulated; 10,000 sweeps (first 20% discarded).
The posterior summary statistics are:

For purposes of comparison, β was also estimated in the usual way (OLS). The estimated specification:

r=(Underscript[1.132, (0.186)]) (Subscript[r, m])+e.

Here are the posteriors:

Extensions

Variation in c

Usual way of analyzing variation in a liquidity proxy:

1. Estimate liquidity proxy over short subsamples.

2. Use the estimates as dependent variables in a subsequent regression.

E.g., Pastor and Stambaugh construct monthly γ estimates where each estimate is based on all the trading days in the month, and then estimate liquidity betas from a portfolio variant of  (Subscript[(OverHat[γ]), i t])=(Subscript[α, i])+(Subscript[β, i])(Subscript[(OverHat[γ]), m t])+(Subscript[e, i t])

Small-sample Gibbs estimates of c don't have nice properties. They are biased. In a small sample, the posterior is dominated by the prior. The prior for c is typically half-normal:

Plot[PDF[NormalDistribution[0,.10],x],{x,0,.3},AxesOrigin->{0,0},PlotRange->All]

ExpectedValue[Abs[x],NormalDistribution[0,.1],x]

It is better to model variation in c directly, incorporating the functional form of c into the Gibbs estimation. Here are some scenarios.

Time-series variation

(Subscript[m, t])=(Subscript[m, t-1])+(Subscript[u, t])
(Subscript[p, t])=(Subscript[m, t])+(Subscript[c, t])(Subscript[q, t]) where (Subscript[c, t])=(Subscript[z, t])γ, i.e., a linear function of observed variables (Subscript[z, t])

Parameter and/or data restrictions may be necessary to ensure (Subscript[c, t])>0.

The Gibbs sampler can be modified to include the simulation of γ as an additional step.

Δ (Subscript[p, t])=(Subscript[c, t])(Subscript[q, t])-(Subscript[c, t-1])(Subscript[q, t-1])+(Subscript[u, t]) = ((Subscript[q, t])(Subscript[z, t])-(Subscript[q, t-1])(Subscript[z, t-1]))γ+(Subscript[u, t])

Given ("conditioning on") the qs and zs, this is a normal regresssion model with coefficient vector γ.
(Subscript[z, t]) might include indicators for event days, etc. ("Are trading costs higher on earnings announcement days?")

Cross-sectional variation

(Subscript[m, i t])=(Subscript[m, i ,t-1])+(Subscript[u, i t])
(Subscript[p, i t])=(Subscript[m, i t])+(Subscript[c, i t])(Subscript[q, ιt]) where (Subscript[c, ιt])=(Subscript[z, i t])(Subscript[γ, i])

Here, z may include cross-sectional variables (e.g., market capitalization) etc.
If the (Subscript[u, i t]) (or (Subscript[q, i t])) are cross-sectionally dependent, the Δ (Subscript[p, t]) regressions should be estimated as a SUR system.
We can avoid this if we impound all cross-sectional dependence in a market factor:

(Subscript[m, i t])=(Subscript[m, i ,t-1])+(Subscript[β, i])(Subscript[M, i t])+(Subscript[u, i t])

Latent variation

(Subscript[p, i t])=(Subscript[m, i t])+(Subscript[c, i t])(Subscript[q, ιt]) where (Subscript[c, i t ])=(Subscript[(OverTilde[z]), t])(Subscript[γ, i])

Here, (Subscript[(OverTilde[z]), t]) is interpreted as an unobserved (latent) liquidity factor that is common to all stocks at given time. We can estimate it because we can write the price change as:

Δ (Subscript[p, i t])=(Subscript[c, i t])(Subscript[q, i t])-(Subscript[c, i t-1])(Subscript[q, i t-1])+(Subscript[u, i t]) =(Subscript[(OverTilde[z]), t])((Subscript[γ, i])(Subscript[q, i t]))-(Subscript[(OverTilde[z]), t-1])((Subscript[γ, i])(Subscript[q, i t-1]))+(Subscript[u, i t])

Written in this fashion, conditional on everything else, (Subscript[(OverTilde[z]), t]) is the coefficient in a cross-sectional regression.

... and more

Bayesian estimation of asset pricing models is well-developed. MCMC estimators (like the Gibbs) are easy to implement because they are modular. Thus, we could estimate an asset pricing model directly on the ms, i.e., returns purged of bid-ask bounce.

Stochastic volatility models are also easily estimated in the Bayesian framework. An SV model could be grafted on the Roll setup, possibly allowing for a link between volatility and the trading cost.

Stochastic duration models may be useful because information intensity is often modeled as being inversely related to intra-trade or intra-event duration.

Missing data models have been used in asset pricing studies to address the problem that we rarely have a full panel of asset returns over time.

A neat Bayesian trick (aside)

Some of the models described above call for panel estimations. Sometimes there are short-cuts.

Suppose that we have a panel regression (Subscript[y, i t])=(Subscript[X, i t])β+(Subscript[e, i t]) where the (Subscript[e, i t]) are not cross-sectionally correlated (E (Subscript[e, i])(Subscript[e, j]^′)=0 for i≠j). We could estimate this as a panel system:

(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {⋮}, {(Subscript[y, N])}}))=(({{(Subscript[X, 1])}, {(Subscript[X, 2])}, {⋮}, {(Subscript[X, N])}}))β+(({{(Subscript[e, 1])}, {(Subscript[e, 2])}, {⋮}, {(Subscript[e, 3])}}))

where (Subscript[y, 1]) is the column vector ((({{(Subscript[y, 1,1]), (Subscript[y, 1,2]), ⋯, (Subscript[y, 1,T])}})))', etc.

Alternatively, suppose that the prior for β is N((Subscript[μ, 0]),(Subscript[Ω, 0])). We proceed sequentially as follows.

Estimate (Subscript[y, 1])=(Subscript[X, 1])β+(Subscript[e, 1]), obtaining a posterior for β, denoted N((Subscript[μ, 1]),(Subscript[Ω, 1])).

Use N((Subscript[μ, 1]),(Subscript[Ω, 1])) as a prior in estimating (Subscript[y, 2])=(Subscript[X, 2])β+(Subscript[e, 2]), obtaining a posterior for β, denoted N((Subscript[μ, 2]),(Subscript[Ω, 2]))

Use N((Subscript[μ, 2]),(Subscript[Ω, 2])) as a prior in estimating (Subscript[y, 3])=(Subscript[X, 23])β+(Subscript[e, 3]), obtaining a posterior for β, denoted N((Subscript[μ, 3]),(Subscript[Ω, 3]))

... and so on. It is sometimes easier to loop in this fashion, rather than lay out and set up the large matrices necessary to run the SUR.

Further note: you can also form a posterior for a single regression by updating observation-by-observation. This is essentially what a Kalman filter does (sequentially, in a time-series context).

Asymmetric information

The Gibbs sampler has been applied to the classic Roll model, in which the increments in the efficient price, (Subscript[u, t]), are assumed to be independent of the trade direction indicators, (Subscript[q, t]).

Could we roughly check this independence by estimating the correlation between the (Subscript[u, t]) and (Subscript[q, t]) (draw by draw)?  Remember that the estimates are coming from the regression

Δ (Subscript[p, t])=c Δ (Subscript[q, t])+(Subscript[u, t])

In a regression, the estimated residuals are uncorrelated with the explanatory variables by construction.

A generalized Roll model that does allow for the trade effects associated with asymmetric information is:

(Subscript[m, t])=(Subscript[m, t-1])+λ (Subscript[q, t])+(Subscript[u, t])
(Subscript[p, t])=(Subscript[m, t])+c (Subscript[q, t])

This implies

Δ (Subscript[p, t])=(λ+c)(Subscript[q, t])+c  (Subscript[q, t-1])+(Subscript[u, t])

The q draws for this model are slightly more involved. See Hasbrouck (2004).

There is a limitation here. The original Roll model holds its form under time aggregation. For example:

(Subscript[p, t])-(Subscript[p, t-2])=c((Subscript[q, t])-(Subscript[q, t-2]))+(Subscript[u, t])+(Subscript[u, t-1])

That is, the bid-ask bounce term is determined by the qs are the endpoints of the observation interval, and the aggregate disturbance has the same form as the original disturbances.

In the generalized Roll model, though:

(Subscript[p, t])-(Subscript[p, t-2])=c((Subscript[q, t])-(Subscript[q, t-2]))+λ((Subscript[q, t])+(Subscript[q, t-1]))+(Subscript[u, t])+(Subscript[u, t-1])

This involves not just the qs at the endpoints, but also the intermediate value (Subscript[q, t-1]). If we were applying the model to closing trade prices, the daily price change between day d-1 and day d would look like:

(Subscript[p, d])-(Subscript[p, d-1])=c((Subscript[q, N])-(Subscript[q, 1]))+λ((Subscript[q, 1])+(Subscript[q, 2])+⋯+(Subscript[q, N]))+((Subscript[u, 1])+(Subscript[u, 2])+⋯+(Subscript[u, N]))

where the subscripts on the q and u run over the N trades on day d. That is, the sum ((Subscript[q, 1])+(Subscript[q, 2])+⋯+(Subscript[q, N])) is the net buy order flow over the day. In principle we could compute the conditional distribution for this sum, and make a Gibbs draw. In practice, though, most of the time we don't know N. This complicates matters enormously.

Volume

Suppose that we're in a Kyle world and (Subscript[y, i]) denotes the net signed volume in the ith market. That is:

Δ (Subscript[p, i])=λ (Subscript[y, i]) where (Subscript[y, i])~N(0,(Subscript[σ, y]^2))

If we had a single observation on Δ p and y, we could estimate λ. It really wouldn't be necessary to know signed volume, though. It would suffice to know the unsigned volume|(Subscript[y, i])|, since |Δ (Subscript[p, i])|=λ | (Subscript[y, i])|.

This intuition, that we can estimate λ from the absolute values, underlies the Amihud illiquidity measure.

Note: in the limit of the Kyle model, where aggregate order flow y is a Wiener process, we might be tempted to interpret |d y| as "volume". This runs into the problem that a Wiener process has unbounded variation.

Aggregation

Things are more complicated under aggregation. Usually we know a day's total trading volume

V=|(Subscript[y, 1])|+|(Subscript[y, 2])|+⋯+|(Subscript[y, N])|

What is f((Subscript[y, 1])+(Subscript[y, 2])+…+(Subscript[y, N])|V)? (And we don't know N.)

Note: If  (Subscript[y, i])~N(0,(Subscript[σ, y]^2)), then |(Subscript[y, i])|/(Subscript[σ, y]) is distributed as a χ variate with one degree of freedom. Note: not a (χ^2) variate. This is unfortunate because (χ^2) variates aggregate nicely, while χ variates do not.

Possible approaches.

Assume that all trades are for 100 shares. Then we're back to inferring (Subscript[q, 1])+(Subscript[q, 2])+…+(Subscript[q, N]). This may be easier.

Easley et al. assume that buys arrive with a Poisson intensity (Subscript[λ, B]) and sells with a Poisson intensity (Subscript[λ, S]). From the aggregation properties of Poisson processes, the trade arrival process has intensity (Subscript[λ, B])+(Subscript[λ, S]). This implies a distribution of volume. Furthermore, the net number of buys is unconditionally distributed as a Skellam variate. What we need, though is the condition distribution of the net number of buys (conditional on the number of trades).
	When the time-aggregation period is sufficiently long that we have a large number of buys and sells, it may be useful to apply the normal approximation to the Poisson distributions. Even when the time period is short, we might still use the normal approximation to generate candidate draws for a Metropolis-Hastings algorithm.

Public (non-trade) information

Δ (Subscript[p, i])=λ (Subscript[y, i])+(Subscript[u, i])=>((| Δ (Subscript[p, i])|)/(| (Subscript[y, i])|))=((| λ (Subscript[y, i])+(Subscript[u, i])|)/(| (Subscript[y, i])|))≤((λ|(Subscript[y, i])|+|(Subscript[u, i])|)/(| (Subscript[y, i])|))=λ+((|(Subscript[u, i])|)/(| (Subscript[y, i])|))

This suggests that I will be a downward-biased estimate of λ.

And finally: a caution

The inferences that we draw from a statistical model depend on the model and the data. We can get stronger inferences from using a tighter model or getting more data (not just more observations, but more kinds of variables).

Sometimes we require data that is either inherently latent or is simply unobtainable for practical reasons. MCMC methods allow us to estimate these data, but only if we are willing to specify a stronger stronger model. That's the trade-off.

The caution is that it is very easy to specify a fanciful model loaded with latent data. If the mixing and convergence properties of the simulation estimator are poor, the model can appear to be well-identified.

Addendum: Transformations of the Gamma distribution

When a=n/2 and b=1/2, the Gamma becomes:

PDF[ChiSquareDistribution[n],x]

We sometimes need to determine the pdf's of (Subscript[σ, u]^2) and/or (Subscript[σ, u]).

If (1/(Subscript[σ, u]^2)) is Gamma, then what is the pdf of (Subscript[σ, u]^2)?

Recall that if y=g(x), then f(y)=f((g^(-1))(y))|(g^(-1′))(y)|

gRule ={g[x_]:>(x^(-1)),gi[y_]:>(y^(-1))};

fy=PDF[GammaDistribution[a,λ],x]*D[-gi[y]/.gRule,y]/.λ->1/b/.x->gi[y]/.gRule//Simplify

Verify that this integrates to unity:

Integrate[fy,{y,0,∞},Assumptions->{a∈Reals, b∈Reals, a>0, b>0}]

Compute the expectation:

Integrate[y fy,{y,0,∞},Assumptions->{a∈Reals, b∈Reals, a>1, b>0}]

If (1/(Subscript[σ, u]^2)) is Gamma, then what is the pdf of (Subscript[σ, u]) (the standard deviation)?

gRuleSD ={g[x_]:>(x^(-1/2)),gi[y_]:>(y^(-2))};

fy=PDF[GammaDistribution[a,λ],x]*D[-gi[y]/.gRuleSD,y]/.λ->1/b/.x->gi[y]/.gRuleSD//Simplify

Verify that this integrates to unity:

Integrate[fy,{y,0,∞},Assumptions->{a∈Reals, b∈Reals, a>0, b>0}]

Compute the expectation:

Integrate[y fy,{y,0,∞},Assumptions->{a∈Reals, b∈Reals, a>1, b>0}]

References

Ball, Clifford A. and Tarun Chordia, 2001, True spreads and equilibrium prices. Journal of Finance 56, 1801-35.
Carlin, Bradley P.  and Thomas A. Louis, 2000,  Bayes and Empirical Bayes Methods for Data Analysis. London, Chapman and Hall.
Geweke, John, 2005,  Contemporary Bayesian Statistics and Econometrics. New York, John Wiley and Sons.
Hasbrouck, Joel, 1999, Security bid/ask dynamics with discreteness and clustering: Simple strategies for modeling and estimation. Journal of Financial Markets 2, 1-28; doi:10.1016/S1386-4181(98)00008-1 
Hasbrouck, Joel, 2004, Liquidity in the futures pits: Inferring market dynamics from incomplete data. Journal of Financial and Quantitative Analysis 39.
Hasbrouck, J. (2009). Trading costs and returns for US equities: estimating effective costs from daily data. Journal of Finance 64(3): 1445-1477 DOI link: http://dx.doi.org/10.1111/j.1540-6261.2009.01469.x 
Kim, Chang-Jin and Charles R. Nelson, 2000,  State-space models with regime switching. Cambridge, Massachusetts, MIT Press.
Lancaster, Tony, 2004,  An Introduction to Modern Bayesian Econometrics. Malden (MA), Blackwell Publishing.
Shephard, Neil, 2005,  Stochastic Volatility. Oxford, Oxford University Press.
Tanner, Martin A., 1996,  Tools for Statistical Inference. New York, Springer-Verlag.