We present an unsupervised, comprehensive methodology for the construction of financial risk models. We o↵ er qualitative comments on incremental functionality and quantitative measures of superior performance of component and mixture dynamic linear models relative to alternative models. We apply our methodology to a high dimensional stream of daily closing prices for approximately 7,000 US traded stocks, ADRs, and ETFs for the most recent 10 years. Our methodology automatically extracts an evolving set of explanatory time series from the data stream; maintains and updates parameter distributions for component dynamic linear models as the explanatory time series evolve; and, ultimately specifies time-varying asset specific mixture models. Our methodology utilizes a hierarchical Bayesian approach for the specification of component model parameter distributions and for the specification of the mixing weights in the final model. Our approach is insensitive to the exact number of factors, and “e↵ ectively” sparse, as irrelevant factors (time series of pure noise) yield posterior parameter distributions with high density around zero. The statistical models obtained serve a variety of purposes, including: outlier detection; portfolio construction; and risk forecasting.
INTRODUCTION We propose a time varying Bayesian statistical model for individual stock returns depicted in Figure 1. Our goal is to accurately model the high dimensional probability distribution underlying stock returns. We demonstrate success by constructing better performing investment portfolios, where the performance goal is to minimize the standard deviation of returns. Investment professionals seek to minimize the variation in investment outcomes because doing so results in higher economic utility for their clients. To illustrate this point, consider which of two pension scenarios with equal expected value is preferred: one that pays 80% salary with probability 1; or, the other that pays 120% salary with probability 1/2 (if things “go well”) and pays 40% salary with probability 1/2 (if things “go poorly”). The economic concept of declining marginal utility, expressed mathematically with concave utility functions U(E(x)) > E(U(x)), implies that given investment scenarios with equal expected return, individuals prefer the scenario with lowest variation. Portfolio managers are concerned with generating acceptable returns with minimal risk; and, risk managers monitor the portfolio managers, verifying financial risks remain within authorized limits. Risk models, defining the n ⇥ n asset covariance matrix ⌃ and precision matrix ⌃ 1, are used by portfolio managers in conjunction with expected return vectors ↵ to construct optimal portfolios weights ⌃ 1↵ ; and, are used by risk managers given portfolio weights w to forecast portfolio variance wT⌃ w.
We describe the construction of a set of Bayesian switching state-space models and qualitatively analyze the on-line behavior of the various component models and mixture models. We focus our discussion on the qualitative aspects then conclude by providing quantitative measures of our model’s superior performance relative to competing models. We look at the tradeo↵ of adaption rate and stability of parameter estimates, evaluating model responsiveness with both synthetic and real data series. We show that dynamic linear models are robust 1 with respect to pure noise explanatory variables, appropriately generating param1We use the term robust to describe a desirable trait where an estimation method is stable in the presence of outlier observations and irrelevant explanatory variables (noise). eter distributions very concentrated around zero. We comment on the behavior of component models relative to the mixture consensus. We illustrate component and mixture response to outlier observations. In periods with ambiguous posterior model probabilities, we describe the di↵ usive impact to the mixture distribution; and, we note surprisingly distinct behavior when an outlier component is selected with near certainty. We find unexpectedly similar updating behavior across a range of component models, bringing into question the necessity of more than two components. We inspect the impact of implementing intervention in the estimation process by greatly inflating the variance of a stock’s posterior distribution subsequent to a merger event. The intervention is shown to result in extremely rapid convergence to new dynamics, while the same model without intervention is shown to maintain bias for an unacceptably long period. Lastly, we compare our two component mixture model against several alternatives. Analyzing results in Table 1, we show the positive impact of regularization provided by the Bayesian framework relative to PCA, and further improvement of the mixture model as compared to the single process model. 2 2.1
BACKGROUND
RISK MODELS High dimensional statistical models, central to modern portfolio management, present significant technical challenge in their construction. There is extensive literature on the topic of constructing factor models or risk models as they are interchangeably known by practitioners. Consider a matrix of observations representing log price returns
2 r1,1 · · · X = 64 ... . . .
rn,1 · · · r1,t 3 .
.. 75 rn,t ri,j = log ✓ pi,j ◆ pi,j 1 , , for assets i 2 1 . . . n, trading days j 2 1 . . . t, and end of day prices pi,j . The general approach to financial risk modeling specifies the covariance of asset returns as a structured matrix. A factor model with p explanatory time series is specified for the n ⇥ t matrix X with an n ⇥ p matrix of common factor loadings L, a p ⇥ t matrix of common factor returns time series F, and an n ⇥ t matrix of residual error time series ✏:
X = LF + ✏ .
An orthogonal factor model
Ft Gt ⇡ ↵ ↵ n,t n,t
Vn,t vn,t
A Yn,t ✓ n,t wn,t W
N
T (1) (2) (3) .
(4)
By construction, risk models based on principal
component and singular value decomposition
COMMON FACTORS
In pursing an unsupervised methodology, we utilize an
SVD based approach to identifying characteristic time
series. SVD identifies common variation, both
rowwise and column-wise, in a matrix of data
We exploit the fact that SVD methodically extracts common variation, ordered by magnitude. The intuition is that “pervasive” sources of return important to modeling portfolio level return characteristics will be reliably captured in the first several factors. Points of concern for some practitioners with regards to the use of SVD or PCA include: 1. the factors are not readily identifiable (to the human analyst); 2. the factors can be permuted in order and sign for data samples adjacent in time; and 3. it’s not clear how many factors to “keep”. With respect to the first concern of (human) identifiability, we reiterate our goal is an unsupervised process yielding a high dimensional statistical model with adequate explanatory power as opposed to semantically meaningful groupings. Examination of assets with significant weight in the characteristic portfolios typically yields meaningful portfolio themes (Johnson and
Wichern, 1998, Ex. 8.5). With respect to the second
concern, the rotation and scaling of factors in di↵
erent sample periods, our application incorporates the
method of
The Bayesian dynamic linear model (DLM)
framework elegantly addresses our need to process streams
of data. DLMs are state space models very similar to
Kalman filters
2 t,n 3 75 , unknown diagonal matrix composed of n-asset specific variance scales at time t; • ⌫ t, n ⇥ 1 vector of unknown observation errors at time t; • ⇥ t = ⇥✓ t,1, . . . , ✓ t,n⇤, a p⇥ n unknown state matrix whose columns are factor loadings of individual assets on common factor returns at time t; and, • ⌦ t = ⇥! t,1, . . . , ! t,n⇤, a p ⇥ n unknown evolution errors matrix applied to the state matrix at time t.
Yt = ⇥ tTFt + ⌫ t, ⌫ t ⇠ N[0, Vt⌃ t] . (6)
tivariate normal, with mean vector 0 and unknown variance Vt⌃ t.
⇥ t = Gt⇥ t 1 + ⌦ t, ⌦ t ⇠ N[0, Wt, ⌃ t] . (7)
⌦ t is matrix normal
evolution of factor loadings across the factors i 2 1 . . . p and across the assets j 2 1 . . . n.
ern, 1998) in Equation 3 to DLM notation of
in Equation 6, the static regression coe cients L now evolve with time in Equation 7 as ⇥ tT.
As typical with a Bayesian approach, our process begins with a prior distribution reflecting our belief about the unknown state parameter matrix ⇥ 0 and the unknown variance scale matrix ⌃ 0 before data arrives. Our initial belief is expressed as a matrix normal/inverse Wishart distribution: (⇥ 0, ⌃ 0) ⇠
NW 01[m0, C0, S0] , (8) with mean matrix m0, left variance matrix C0, right variance matrix S0, and degrees of freedom 0. We allow our estimate of observational variance to vary over time by decaying our sample variance degrees of freedom parameter t 1 immediately before computation of the prior distribution Equation 10.
matrix T distribution: ⇥ 0 ⇠
T 0 [m0, C0, S0] . (9) Let Dt = ⇥Yt . . . Y0⇤ refer to the information available subsequent to observing Yt. The conjugate parameter distributions are updated as follows.
(⇥ t, ⌃ t|Dt 1) ⇠
NW t1 1 [at, Rt, St 1] , (10) where at = Gtmt 1 and Rt = GtCt 1GtT + Wt.
vector Ft: (Yt|Dt 1) ⇠ T t 1 [ft, QtSt 1] , (11) where ft = atTFt and Qt = Vt + FtTRtFt .
is observed. Therefore, we accommodate
random regression vectors Ft
known is (Yt|Dt 1) ⇠ T t 1 [fˆt, QˆtSt 1], where the moment parameters of the multivariate T forecast distribution are now ˆft = atTµFt and QˆtSt 1 = T
Vt + µTFt RtµFt + tr (Rt⌃ Ft ) St 1 + at ⌃ Ft at.
(⇥ t, ⌃ t|Dt) ⇠
NW t1[mt, Ct, St] , (12) with mt = at + AtetT, Ct = Rt AtAtTQt, t = t 1 + 1 and St = t 1 ⇥ t 1St 1 + etetT/Qt⇤ where At = RtFt/Qt and et = Yt ft. 2.4
UNIVARIATE DLMS
DLMs. Setting the number of assets n = 1, results
for univariate DLMs immediately follow.
With class II models, there are |A|t potential model histories for each asset. We avoid this explosion in model sequences by considering only two periods, t 1 and t, thereby limiting distinct model sequences in our mixtures to |A|2. As each asset has its own history, no longer sharing common scale free posterior variance Ct, our mixture models are asset specific, forcing the use of univariate component DLMs. Parameters 1 ⇥ 1 in univariate DLMs are now displayed with scalar notation. To avoid clutter, we omit implicit asset subscripts.
Inference with multi-process models is based upon manipulation of various model probabilities: the posterior model probabilities for the last model pt 1(↵ t 1); the prior model probabilities for the current model ⇡ (↵ t); and, the model sequence likelihoods p(Yt|↵ t, ↵ t 1, Dt 1). Posterior model sequence probabilities for current model ↵ t and last model ↵ t 1 upon observing Yt are: pt(↵ t, ↵ t 1) = Pr[↵ t, ↵ t 1|Dt] / pt 1(↵ t 1)⇡ (↵ t)p(Yt|↵ t, ↵ t 1, Dt 1) . (13) The unconditional posterior parameter distributions are computed as mixtures of the |A|2 component DLM sequences
k p(✓ t|Dt) = X k
X ↵ t=1 ↵ t 1=1 The posterior model probabilities are
pt(✓ t|↵ t, ↵ t 1, Dt)pt(↵ t, ↵ t 1) . pt(↵ t) = Pr[↵ t|Dt] = pt(↵ t, ↵ t 1) . k
X The posterior probabilities for the last model ↵ t 1 given the current model ↵ t and information Dt are (14) (15) .
(16)
After each time step, the posterior mixture
distribution for each component DLM is approximated with
an analytic distribution using the methodology
described in
St(↵ t) 1 =
1 Xk pt(↵ t, ↵ t 1) pt(↵ t) ↵ t 1=1 St(↵ t, ↵ t 1) .
(17) The weights for computing the moments of the KL divergence minimizing approximation to the posterior distribution are: p⇤t (↵ t 1) =
St(↵ t) pt(↵ t, ↵ t 1) pt(↵ t) St(↵ t, ↵ t 1) .
(18) The mean vector mt(↵ t) for DLM ↵ t is: mt(↵ t) = p⇤t (↵ t 1)mt(↵ t, ↵ t 1) . (19) k
X The variance matrix Ct(↵ t) for DLM ↵ t is: Ct(↵ t) = k
X ↵ t 1=1
p⇤t (↵ t 1) n Ct(↵ t, ↵ t 1) + [mt(↵ t) [mt(↵ t) mt(↵ t, ↵ t 1)] ⇥
T o
mt(↵ t, ↵ t 1)]
. (20)
The degrees of freedom parameter, t(↵ t) is important
to the KL minimization. Intuitively, if the component
DLMs are in discord, the resulting mixture may be
described as “fat-tailed”, and the precision of the
unknown variance scale parameter reduced. We compute
t(↵ t) using the procedure described by
Q2
Q3
Q4 2012 Q1
Q2
Q3
Q4 (a) Synthetic scenario 1: abrupt increase in factor loading ✓ t, units are percent per annum. Observations Yt are synthesized returns using SPY (S&P 500 ETF) until March 30, 2012; and, 2⇥ SPY thereafter. Note “true” factor loading doubles from approximately 10 to 20. (b) Synthetic scenario 2: abrupt decrease in factor loading ✓ t, units are percent per annum. Observations Yt are synthesized returns using SPY until March 30, 2012; and, AGG (Barclays Aggregate Bond ETF) thereafter. Note “true” factor loading drops from approximately 10 to 0.
APPLICATION DESIGN
Our application enables a proprietary trading group at a financial services firm to better understand the aggregate behavior of stock portfolios. The models provide a statistical framework required for constructing portfolios, assessing portfolio risk, and clustering assets. The assets of primary interest are the common shares of the largest 1000 - 2000 companies in the US stock market. The group focuses on larger companies because the liquidity of larger stocks generally makes them cheaper and easier to transact. The group’s strategies include short selling: borrowing stock, selling borrowed shares; and, attempting to profit by repurchasing the shares at a lower price. Larger stocks are generally easier to borrow. 3.2
In implementing our application, one of the first issues
encountered is a machine learning classic, the
biasvariance tradeo↵ [Ch. 2.9]
We identify two universes of assets: a relatively narrow set that will be used to construct explanatory time series; and, an all inclusive set for which we will generate factor loading and residual volatility estimates.
It is desirable that the assets used to construct the
common factor returns trade frequently and with
adequate liquidity to minimize pricing errors. We also
avoid survivor bias, the methodological error of
omitting companies no longer in existence at the time a
historical analysis is performed. We eliminate this
hazard by defining our common factor estimation universe
using daily exchange traded fund (ETF) create /
redeem portfolios for the last 10 years. Brokers create
ETF shares (in units of 50,000 ETF shares) by
delivering a defined portfolio of stock and cash in exchange
for ETF shares; or, they redeem ETF shares and
receive the defined portfolio of stock and cash. Given the
create / redeem definitions that were used as the basis
for large transactions during the historical period, the
assets in an ETF portfolio represent an institutionally
held, survivor bias free, tradeable universe. Its likely
the component shares were readily available to borrow,
as the component shares were held by custodian banks
for the ETF shareholders. Our data base permits us
to obtain data for surviving and extinct stocks. The
ETF we select for our factor estimation universe is the
Vanguard Total Stock Market ETF (VTI)
DATA PREPARATION
Data preparation involves constructing the artifacts
demanded by § 2.3: Yt, Ft, and Gt. Each trading day,
using price, dividend, and corporate action data for all
7000 - 8000 stocks, ADRs, and ETFs in our US pricing
data base
1
ical returns matrix rt = ⇥Yt T +1 . . . Yt⇤ ⌃ ˆ t 2 where
⌃ ˆ t is the diagonal matrix of sample variance for the
period t T + 1 to T . Using
1 adjacent data windows, rt 1 = ⇥Yt T . . . Yt 1⇤ ⌃ ˆ t 21
1
and rt = ⇥Yt T +1 . . . Yt⇤ ⌃ ˆ t 2 are then used to
compute Gt as described in
Gt = FtFtT 1
FtFtT 1 .
(21) We are a little flexible with the notation in Equation 21, where Ft and Ft 1 are p⇥ (T 1) sub-matrices representing time aligned subsets of two factor return matrices, the scaled right singular vectors obtained from the decomposition of rt 1 and rt. Elsewhere, viz. Equation 11, Ft refers to a p ⇥ 1 dynamic regression vector, the right most column of the transposed and scaled right singular vectors, corresponding to the desired vector of common factor returns for day t. 3.5
MODEL COMPONENTS
The component DLMs in our mixture model share
observations Yt, common factor scores Ft, and state
evolution matrices Gt. The component DLMs are di↵
erentiated by the variance parameters: the observational
variance scale Vt, and the evolution variance matrix
Wt. We construct a set of component DLMs following
the approach of
For component DLMs that will accommodate “typical” observation variance, we set Vt = 1; for component DLMs that will accommodate outlier observations, we set Vt = 100. For the evolution variance, we similarly select a base rate of evolution Wt = .00001; and, inflate Wt by a factor of 100 and 1002 to permit increasingly rapid changes in the factor loadings. 3.6
OUTPUT The format of the output risk model will be a n ⇥ p factor loading matrix and a n ⇥ 1 residual volatility vector. These p + 1 numeric attributes for n stocks are stored in various formats for subsequent use throughout the organization. 4
EVALUATION
Factor 1 Factor 2 Factor 3 Factor 4 Factor 5
N(0,1) noise 600 500 400 300 200 100 0 -0.010 -0.005 =12 .00001
(c) 10 uWre
W Aug
Mixture 4 V=1 W=.00001
V=1 W=.001 =0 .00001 2 V=1 W=.1
V=100 W=.00001 Jul Aug Sep Oct
2013 W160 150
=.1 =W200.
0=0.010001
(e)
Sep =W180.0=0.000011 190 IBM 170
SPY
2013
W=.00001
Jul Aug Sep Oct Nov Dec
2013 (g) Price histories for IBM and SPY, units are USD.
Sep Oct Nov Figure 4: Multi-Process Models. See § 4.3 for discussion. Figure 4(a),(c),(e) units are percent per annum. (h) Posterior component probabilities (left), color key in (c); mixture model degrees of freedom t (right), white line.
Sep 2013 Jul
Aug
Oct
Nov Dec Aug Nov
Dec Nov
Dec (b) (d) (f) 1.0 0.8
Oct 0.6 0.4 0.2 0.0
Nov
12
10
8
6
4
2
regression parameter distributions densely
surrounding zero given su cient observations. Further, a
linear combination of independent DLMs is a DLM
In our application, we focus on aggregate (portfolio
level) forecasts, therefore it is extremely important to
identify common sources of variation that may appear
insignificant at the individual asset level.
MIXTURE DYNAMICS Figure 4 shows the price movement and model response for IBM during the second half of 2013. In Figure 4(g), the price histories for IBM and the S&P 500 ETF SPY are displayed. Note the sharp drop in IBM’s price on October 17, 2013 corresponding to an earnings announcement. This event is helpful in understanding the interaction of components in our multi-process models. We construct three multi-process models, the simplest of which is a two component mixture (the “base model”), comprised of a standard component DLM to handle the majority of the observations Yt, and an outlier component DLM. We specify common evolution variance Wt = .00001 I; observation variance Vt = 1 for the standard component; and observation variance Vt = 100 for the outlier component. The base model and component estimates for the first factor loading m1,t appear in Figure 4(a) and magnified in Figure 4(b). We specify a three component mixture (the “adaptive model”) by adding a component DLM with inflated evolution variance Wt = .001 I. The adaptive model and component estimates for the first factor loading m1,t appear in Figure 4(c) and magni2A basis point (bp) is 10, 000 1. fied in Figure 4(d). Finally, we specify a four component mixture (the “very adaptive model”) by adding a component with further inflated evolution variance Wt = .1 I. The very adaptive model and component estimates for the first factor loading m1,t appear in Figure 4(e) and magnified in Figure 4(f). The posterior component model probabilities for the very adaptive model appear as a bar chart in Figure 4(h), where the bottom bar corresponds to the probability of the outlier component, the second bar corresponds to the very adaptive component DLM, the third bar corresponds to the adaptive component DLM, and the top bar corresponds to the base component DLM. We specified the fixed DLM selection (prior) probabilities as .01543, .00887, .0887, .887 for the components {outlier, very adaptive, adaptive, base} respectively.
Note several occurrences where the posterior probability of an outlier observation significantly exceeds the DLM selection probability. The white line in Figure 4(h) corresponds to the degrees of freedom parameter t for the T-distribution that approximates the mixture model’s posterior parameter distribution. 4.3
QUALITATIVE COMMENTS To supplement the more analytically precise discussion in § 2.5, we make the following qualitative comments as to the interaction of the mixture components: • the mixture in Figure 4(e) with larger evolution variance adapts faster; the mixture in Figure 4(a) with smaller evolution variance appears smoother; • time t component posteriors are 1-period departures from the t 1 consensus, see Figure 4(a), (b), (c), (d), and (e); • an outlier component “ignores” current observations Yt and forecast error |et|, responding to the t 1 posterior consensus mt 1, see Figure 4(b) and (d); • in periods of noise, the other components return to the outlier component’s estimate with 1-period lag, see left-hand side of Figure 4(b) and (d); • in periods of level change, the outlier follows the other components’ estimate with 1-period lag, see right-hand side of Figure 4(b) and (d); • when the posterior probability of the outlier component spikes up, the degrees of freedom parameter t usually drops, reducing the precision of the variance scale estimate St, see Figure 4(h); • however, when the outlier component is selected with very high probability, there is no impact to t as the observation is ignored, see October 17, 2 dlm mix + intervention Q2 2013 in Figure 4(h), noting that the white line does not drop when Pr{ outlier } ⇡ 1; • except for the very adaptive component, the response does not vary with Wt, see Figure 4(d) and (f), where W=.001 and W=.00001 responses are nearly identical. 4.4
UNRESPONSIVENESS TO Wt The phenomenon we find most surprising is the insensitivity of the components to Wt below a certain threshold. Digging further into this phenomenon, in a mixture, the various component models view of the t 1 posterior parameter distribution are very similar. Thinking about univariate DLMs, and assuming for discussion F = 1 and G = 1, the magnitude of the adaptive scalar At = Rt/Qt = (Ct 1 + Wt) / (Ct 1 + Wt + Vt). When Wt ⌧ Ct 1, as describes our situation, At ⇡ Ct 1/ (Ct 1 + Vt) as seen in Figure 4(d) and (f). Only when Wt is significant relative to Ct 1 does response vary noticeably. 4.5
INTERVENTION
Our discussion of IBM focused on the mixture
models’ processing of unusual observations. We now
explore an example where the data generating process
changes abruptly, similar to our synthetic illustrations
in Figure 2. In April 2013, the acquisition of Life
Technologies Corporation by Thermo Fisher Scientific
Inc. was announced. The stock’s sensitivity to the
market, as expressed in its first common factor
loading, dropped abruptly, as shown in Figure 5. While
our goal is an unsupervised estimation processes, the
Bayesian DLM framework facilitates structured
intervention when necessary. For one of the mixture models
in Figure 5b, we intervene and inflate the prior
parameter variance following the April 15th announcement,
Rt++ = Gt (Ct 1 + I) GT + Wt, where the identity
matrix reflects the increased uncertainty in the
parameter distribution relative to the usual prior variance in
Equation 10. When subsequent updates occur, the
DLM with the inflated prior variance adapts to the
new dynamics rapidly and satisfactorily.
To access the performance of our two component
mixture model, we construct daily portfolios from the VTI
universe for the most recent ten years, May 2004 to
April 2014. The number of trading days during this
period was 2,516. In Table 1, we report realized
outof-sample volatility and the standard error (se) of the
volatility measure for two strategies: global minimum
variance (GMV); and, maximum Sharpe ratio (MSR)
CONCLUSION The ability to integrate SVD with Bayesian methods allows our application to process large data streams in an unsupervised fashion. We demonstrate that a two component multi-process model achieved better reduction in volatility than alternative models. The two component model out-performed alternative models including a single process model. We find the robustness of Bayesian DLMs with respect to noise inputs of great practical value, allowing us to favor inclusion of factors, potentially capturing pervasive sources of common movement important to aggregate forecasting. The inclusion of an outlier model adds great functionality, delivering robustness to the estimation process. The insensitivity of the mixture models to multiple evolution variance values leads us to favor mixtures of just two components, a typical evolution variance value for both components, and an inflated observation variance in the outlier component. We would recommend generating several models of varying adaptiveness, evaluating the variance-bias tradeo↵ in light of a user’s specific situation. We favor the judicious use of intervention for events such as mergers.
We would like to explore using news feeds to systematically intervene for events known to impact an assets dynamics.
Applied
G. Connor and R.A. Korajczyk. Risk and return in an