Dynamic Portfolio Optimization

A DISSERTATION
SUBMITTED TO THE FACULTY OF THE UNIVERSITY OF
MINNESOTA BY

Lei Yuan

IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE

DEGREE OF DOCTOR OF PHILOSOPHY

Fadil Santosa, Daniel Mitchell

June, 2021
© Lei Yuan 2021
Acknowledgements

I would like to thank my advisor Fadil Santosa for accepting me as his student, guiding
me through the graduate program, and helping me start my professional career. Working
with him has always been an encouraging experience. I would also like to thank my
co-advisor Daniel Mitchell, who has provided me with valuable research insights and
advice. I am grateful to Gary Hatfield and Securian Financial for the support to my
research and the opportunity to gain industry experience. Finally, I would like to thank
my family and friends, without whom I could not have made it this far.

i
Dedication

To my wife and parents.

ii
 Abstract

The problem of portfolio optimization has attracted broad attention from the academia
and financial industry. In this thesis, we will explore a variety of mathematical models
related to portfolio optimization. We will derive equations that describe price dynamics
of risky assets, and study portfolio optimization problems based on these equations.

iii
Contents

Acknowledgements i

Dedication ii

Abstract iii

List of Figures vi

1 Introduction 1

2 Robust Portfolio Optimization 4

2.1 Robust Portfolio Selection Model . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Linear Regression and Uncertainty Sets . . . . . . . . . . . . . . . . . . 6
2.3 Joint Ellipsoidal Uncertainty Set . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Numerical Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Bayesian Methods 13
3.1 The Black-Litterman Model . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Extension of Black-Litterman Model . . . . . . . . . . . . . . . . . . . . 16
3.3 Robust Bayesian Allocation . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4 Numerical Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4 Dynamic Portfolio Optimization 22

4.1 Robust Dynamic Programming . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Merton’s portfolio problem . . . . . . . . . . . . . . . . . . . . . . . . . 26

iv
5 Dynamic Portfolio Optimization with Bayesian Update - Discrete Model 31
5.1 A Discrete Model with an Uncertain Parameter . . . . . . . . . . . . . . 31
5.2 Transformation to a Dynamic Programming . . . . . . . . . . . . . . . . 34
5.3 Numerical Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

6 Dynamic Portfolio Optimization with Bayesian Update - Continuous

Model 39
6.1 Conjugate Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
6.2 Stochastic Market Model . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.3 Stochastic Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.4 Stochastic Control Problem - Unconstrained Action . . . . . . . . . . . 45
6.5 Stochastic Control Problem - Single Action . . . . . . . . . . . . . . . . 53
6.6 Stochastic Control Problem - Finite Action Set . . . . . . . . . . . . . . 57
6.7 Stochastic Control Problem - Constrained Action . . . . . . . . . . . . . 74
6.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

7 Conclusion 87

References 88

Appendix A. 92
A.1 Derivation of HJB Equation . . . . . . . . . . . . . . . . . . . . . . . . . 92

v
List of Figures

1 Summary statistics of asset returns . . . . . . . . . . . . . . . . . . . . . 10

2 Portfolio allocations using GI and Lu models . . . . . . . . . . . . . . . 11
3 Summary statistics and portfolio allocation using RBA model . . . . . . 21
4 Optimal action tree (no short sell allowed) . . . . . . . . . . . . . . . . . 37
5 Optimal action tree (short sell allowed) . . . . . . . . . . . . . . . . . . 38
6 Numerical solution of A at some specific times . . . . . . . . . . . . . . 46
7 Numerical solution of B(t, x) at some specific times . . . . . . . . . . . . 47
8 Numerical derivative of B(t, x) w.r.t x at some specific times . . . . . . 48
9 Numerical solution of B(t, x) at some specific times (α = 0.5) . . . . . . 54
10 Numerical solution of B(t, x) and optimal control . . . . . . . . . . . . . 61
11 Comparison of the numerical solution with analytical solutions to single
action problems (A = {0} and A = {1}) . . . . . . . . . . . . . . . . . . 62
12 Numerical derivatives of B(t, x) . . . . . . . . . . . . . . . . . . . . . . . 64
13 Numerical solution of B(t, x) and optimal control . . . . . . . . . . . . . 66
14 Comparison of the numerical solution with analytical solutions to single
action problems (A = {−1} and A = {1}) . . . . . . . . . . . . . . . . . 67
15 Numerical derivatives of B(t, x) . . . . . . . . . . . . . . . . . . . . . . . 69
16 Numerical solution of B(t, x) and optimal control . . . . . . . . . . . . . 71
17 Comparison of the numerical solution with analytical solutions to single
action problems (A = {−1} and A = {1}) . . . . . . . . . . . . . . . . . 72
18 Numerical derivatives of B(t, x) . . . . . . . . . . . . . . . . . . . . . . . 73
19 Numerical solution of B(t, x) and optimal control . . . . . . . . . . . . . 76
20 Comparison of the numerical solution with analytical solutions to single
action problems (A = {0} and A = {1}) . . . . . . . . . . . . . . . . . . 77

vi
21 Numerical derivatives of B(t, x) . . . . . . . . . . . . . . . . . . . . . . . 78
22 Numerical solution of B(t, x) and optimal control . . . . . . . . . . . . . 80
23 Comparison of the numerical solution with analytical solutions to single
action problems (A = {−1} and A = {1}) . . . . . . . . . . . . . . . . . 81
24 Numerical derivatives of B(t, x) . . . . . . . . . . . . . . . . . . . . . . . 82
25 Evolution of optimal control in time (A = [−1, 1]) . . . . . . . . . . . . 84
26 Evolution of optimal control in time (A = [−30, 30]) . . . . . . . . . . . 85

vii
Chapter 1

Introduction

Portfolio optimization is the problem of finding the optimal allocation of funds over
a number of available assets in order to maximize the return on the investment while
minimizing the risk. Solutions to this problem is important for both institutions and
individuals. Institutions such as mutual funds and hedge funds make portfolio decision
on a regular basis. Individuals planning for retirement are also concerned about the
implications of their choices. Mathematical models have been developed to address var-
ious issues of relevance for the parties involved.

Makowitz [27] introduced the first mathematical model for portfolio optimization in
the 1950s. Let the random vector r = (r1 , ..., rn )T ∈ Rn denote returns of n assets in the
market, and φ = (φ1 , ..., φn )T ∈ Rn denote the proportions of wealth to be invested in
the n risky assets such that 1T φ = 1. Suppose r has mean µ and covariance matrix Σ.
According to Markowitz’s theory, portfolio optimization is to find the optimal trade-off
between “return” and “risk”, which are represented by the expected return φT µ and
the variance φT Σφ, respectively. Mathematically speaking, given a minimum expected
return µ0 ∈ R, one has the optimization problem

minimize
n
φT Σφ
φ∈R

subject to φT µ ≥ µ0 , (1)

1T φ = 1.

1
 2

The formulation suggests that, given a lower bound on the expected return the in-
vestor is willing to accept, the optimal portfolio can be obtained by minimizing the vari-
ance. If parameters µ and Σ are known, the solution to (1) can be found analytically[29].
It can also be solved numerically under additional constraints, such as no short-selling
and position limits in any single asset.

For any µ0 such that (1) is feasible, the optimal solution φ∗ gives the portfolio al-
location among n risky assets which minimizes variance. In the return-variance space,
the relationship between constraint µ0 and the optimal value of objective function,
φ∗ T Σφ∗ , can be represented by a curve known as the efficient frontier. Intuitively, an
investor should always allocate the capital to stay on the efficient frontier. Based on
Markowitz’s work, Treynor[36, 37], Sharpe[34], and Lintner[25] developed the Capital
Asset Pricing Model (CAPM) for pricing an individual security or portfolio. CAPM
has a profound impact in the academia and practitioners. It has been generalized into a
family of factor-based investing models such as the Fama-French models[11, 12], which
are widely used in today’s financial industry.

The first part of this dissertation addresses the issue of uncertainty in Markowitz
model. Since in real-world applications, parameters such as µ and Σ are not known
exactly and are usually estimated using historical data of financial securities and statis-
tical regression technique, which introduce noise and error to estimated values. Because
the quadratic optimization problem in Markowitz model is sensitive to small perturba-
tion of µ and σ, the result could be problematic[4] and many practitioners turn away
from it. To address this issue, several approaches have been introduced. Goldfarb
and Iyengar[16] developed a robust optimization framework to take possible ranges of
parameters into consideration. Black and Litterman[7, 6] introduced a Bayesian frame-
work to modify statistical estimations with investors’ subjective view. Meucci[31] latter
unified two approaches into a robust Bayesian model.

The second part of the dissertation studies portfolio optimization problems from a
dynamic perspective. As the financial market is ever evolving, investors need to adjust
 3
portfolio allocations according new data and opinion. Thus, portfolio selection can be
viewed as dynamics optimization problems in discrete time or continuous time. Dis-
crete portfolio optimization requires less mathematical foundation and can usually be
solved using discrete dynamic programming. Continuous portfolio optimization, which
was first studied by Merton[28, 30] in the 1970s, sits on stochastic calculus and partial
differential equation. We will review both types of problems and focus on a problem
which is in the intersection of stochastic control and Bayesian theory.

The rest of the dissertation is organized as follow. In chapter 2, we review robust

portfolio optimization models developed to address drawbacks of Markowitz’s model. In
chapter 3, we look at optimization models that borrow ideas from on Bayesian statistics.
In chapter 4, we go back an review the classic continuous-time portfolio optimization
problem proposed by Merton. In chapter 5 and 6, we introduce and study dynamic
portfolio optimization models in discrete-time and in continuous-time.
Chapter 2

Robust Portfolio Optimization

Although Markowitz’s mean-variance model provides a convenient and useful theoretical

framework for portfolio optimization, it has several drawbacks. Broadie[9] has shown
that the optimal portfolio is sensitive to perturbations of market variables µ and Σ.
Since µ and Σ are subject to statistical estimation error, the result of its following op-
timization problem is not reliable. In addition, the model can and often lead to poorly
diversified portfolio by allocating capital to only a few assets. To reduce the sensitivity
of the Markowitz’s optimal portfolios to input error, Goldfarb and Iyengar[16] proposed
portfolio selection models that are robust to parameter uncertainty and estimation er-
rors, and can be efficiently solved as convex optimization problems.

2.1 Robust Portfolio Selection Model

Assume that the market opens for trading at discrete time and has n traded assets. The
vector of asset returns over a single market period is denoted by a column vector r ∈ Rn .
The returns on the assets in different market periods are assumed to be independent.
The single period asset return r ∈ Rn is assumed to be a random variable driven by the
factor model
r = µ + VT f + , (2)

where µ ∈ Rn is the mean return, f ∼ N (0, F) ∈ Rm denotes the normally distributed

return of factors that drives the market with mean 0 and covariance F, V ∈ Rm×n is the

4
 5
deterministic factor loading matrix, and  ∼ N (0, D) ∈ Rn is a normally distributed
residual vector with mean 0 and covariance D. In addition, we assume  is independent
of f , F is positive definite, and D is a positive semi-definite diagonal matrix with
diagonal elements d ∈ Rn (i.e. D = diag(d)). According to the factor model, the asset
returns is also normally distributed with r ∼ N (µ, VT FV + D). Denote the position of
investment as φ ∈ Rn , where the ith component φi represents the fraction of total wealth
invested in the ith asset. The return of the portfolio φ is then normally distributed with
distribution
rφ ∼ N (φT µ, φT (VT FV + D)φ). (3)

In contrast to Markowitz’s model where mean return and the covariance matrix are
constant, Goladfarb and Iyengar (GI) introduced uncertainty sets for D, V, and µ in
order to deal with statistical errors. In GI model, each diagonal element of the covariance
matrix D can take any value in an closed interval and the uncertainty set Sd for D is

Sd = {D : D = diag(d), di ∈ [di , di ] for i = 1, ..., n}. (4)

Each element of mean return vector µ can take values from an interval centered at the
ith element of a constant vector µ0 ∈ Rn

Sm = {µ : µ = µ0 + ξ, |ξi | < γi for i = 1, ..., n}, (5)

where γi ’s are constant real numbers. The uncertainty set of loading matrix V is given
by restricting each column to an elliptic ball w.r.t. a symmetric positive definite matrix
G ∈ Rm×m , i.e.

Sv = {V : V = V0 + W, ||Wi ||g ≤ ρi for i = 1, ..., n}}, (6)

√
where Wi denotes the ith column of W, and ||w||g = wT Gw for any w ∈ Rm . The
choices of uncertainty sets is motivated by the fact that parameters where uncertainty
sets are established around are estimated by linear regression between historical time
series of r and f . Other parameters such as di , di , µ0 , and V0 can be obtained by
multivariate linear regression and by specifying confidence levels.
 6
Combining (3) and uncertainty sets, GI proposed a robust version of Markowitz’s
portfolio optimization problem

minimize
n
max φT (VT FV + D)φ
φ∈R {V∈Sv ,D∈Sd }

s.t. min φT µ ≥ α (7)

{µ∈Sm }

1T φ = 1,
where α ∈ R stands for the worst-cast return level. The objective of the robust min-
imum variance portfolio selection problem (7) is to minimize the worst-case variance
of the portfolio subject to the constraint that the worst-case expected return on the
portfolio is at least α. Since it takes into account all possible values of V, D, and
µ in uncertainty sets, it is less sensitive to parameter perturbation than the original
formulation of Markowitz.

With uncertainty sets Sd , Sm , and Sv described in (4), (5), and (6), the robust
optimization problem can be converted into a Second Order Cone Programming (SOCP)
problem and solved by numerical convex optimization software such as CVX[17].

2.2 Linear Regression and Uncertainty Sets

GI provided a way of getting statistical estimations of uncertainty sets Sd , Sm , and Sv
via multivariate linear regression. Suppose the market data consists of asset returns
{rt : t = 1, ..., p} for p ∈ N periods, and the corresponding factor returns {f t : t =
1, ..., p}. Then the factor model r = µ + VT f +  implies that
m
X
rit = µi + Vji fj + ti (8)
j=1

for 1 ≤ i ≤ n and 1 ≤ t ≤ p. In linear regression analysis, it is standard to assume

the vector of residual returns t ’s are independent and identically distributed normal
random variables, i.e. ti ∼ N (0, σi2 ) for t = 1, 2, ..., p. Let S = [r1 , ..., rp ] ∈ Rn×p be the
matrix of observed asset returns and B = [f 1 , ..., f p ] ∈ Rm×p be the matrix of observed
factor returns. The factor model implies the following linear model for each asset

yi = Axi + i , (9)
 7
where yi = [ri1 , ri2 , ..., rip ]T , A = [1, BT ], xi = [µi , V1i , V2i , ..., Vmi ]T , and i = [e1i , e2i , ..., epi ]
is the residual returns corresponding to asset i. If we assume AT A is invertible, which
is usually the case when p is large enough, the least squares estimation xi of xi is

xi = (AT A)−1 AT yi . (10)

Substituting (9) into (10) yields

xi − xi = (AT A)−1 AT i ∼ N (0, σi (AT A)−1 ). (11)

1
As a result Xi = σi2
(xi − xi )T (AT A)(xi − xi ) is a χ2 random variable with m + 1 degrees
of freedom. Since the true value of σi2 is unknown, we can replace σi2 by its unbiased
||yi −Axi ||22
estimation s2i = p−m−1 and get the resulting random variable
1
Yi = (xi − xi )T (AT A)(xi − xi ), (12)
(m + 1)s2i
which is distributed according to the F −distribution with m + 1 degrees of freedom[1]
in the numerator and p − m − 1 degrees of freedom in the denominator. Now for any
number ω ∈ (0, 1), there is a ω-critical value c(ω) such that the probability of the event
Yi ≤ c(ω) is ω. Equivalently, the set Si (ω) defined below is a ω−confidence set for xi

Si (ω) = {xi : (xi − xi )T (AT A)(xi − xi ) ≤ (m + 1)c(ω)s2i }. (13)

By independent assumption, S(ω) = S1 (ω) × S2 (ω) × ... × Sn (ω) is a ω n -confidence set

for (µ, V). Based on S(ω), we can then get a natural description of an ω n -confidence
uncertainty set for µ by projecting S(ω) along the direction of µ:

Sm (ω) = {µ : µ = µ0 + v, |vi | ≤ γi , i = 1, 2, ..., n}, (14)

q
where µ0 is sample average of historical returns, and γi = (AT A)−1 2
11 c(ω)si . Similarly,
the uncertainty set for the loading matrix is:

Sv (ω) = {V : V = V0 + W, ||Wi ||g ≤ ρi , i = 1, 2, ..., n}, (15)

where V0 comes from the least-square estimation in (10), ||Wi ||g is the elliptic norm
of the ith column of W with respect to the matrix G = BBT − p1 (B1)(B1)T , and
q
ρi = mc(ω)s2i . GI didn’t give a specific description of uncertainty set Sd . It is
understood to be more flexible than Sm and Sv and can be treated as a known parameter
and replaced by empirical estimations.
 8
2.3 Joint Ellipsoidal Uncertainty Set
There are several criticisms of GI on the uncertainty sets. One of them is about the size
of uncertain sets Sv and Sv . Because they are constructed by projecting S(ω) along
the directions of µ and V. The inner maximization problem in (7) optimizes over the
Cartesian product Sm × Sv , which can be substantially larger than the original uncer-
tainty set S(ω). In addition, the actual probability of Sm × Sv is unknown, and the
optimal portfolio tends to be too conservative and poorly diversified.

Lu[26] introduced a joint ellipsoidal uncertainty set to address this issue within GI’s
framework. He also gave a natural description of the joint ellipsoidal uncertainty set
based on multivariate linear regression and an equivalent formulation of problem (7).
Instead of looking at each individual random variable Yi defined by equation (12), Lu
Pn
considered their sum Y = i=1 Yi and defined the joint ellipsoidal uncertainty set
Sm,v (ω) as
n
X (xi − xi )T (AT A)(xi − xi )
Sm,v (ω) = {(µ, V) : Y = ≤ c̃(ω)}, (16)
i=1
(m + 1)s2i

where c̃(ω) is a number depends on ω and the distribution of Y such that the probability
of Y ≤ c̃(ω) is ω. Although this new formulation avoids projecting sets and taking
Cartesian product, one problem that comes up is the the transformed random variable
Y does not have an analytic cumulative distribution function (CDF). As a result, finding
c̃(ω) given ω is not as easy as in GI model, where one can simply look up c(ω) from the
CDF of an appropriate F −distribution. To handle this problem, Lu suggested that one
can either find c̃(ω) by running numerical simulations if the number of assets n is not
large, and approximate the CDF of Y via the Central Limit Theorem if n is large. The
corresponding robust optimization formulation using ellipsoidal uncertainty set is

minimize
n
max φT (VT FV + D)φ
φ∈R {(µ,V)∈Sm,v }

s.t. min φT µ ≥ α, (17)

{(µ,V)∈Sm,v }

1T φ = 1.
 9
Using duality argument, Lu showed the above optimization problem has the same solu-
tion to the following problem for some positive constant θ:

maximize
n
max φT µ − θφT (VT FV + D)φ
φ∈R {(µ,V)∈Sm,v }
(18)
T
s.t. 1 φ = 1.

Lu further proved that problem (18) can be efficiently solved after being transformed
into a SOCP problem and claimed the resulting optimal portfolios are more diversified
than the solutions based on GI model.

2.4 Numerical Study

We present some numerical study result to see how GI and Lu models perform with
real financial data. In the study, we use historical price data of 13 assets to estimate
their mean return vector, covariance matrix, loading matrix, and other model param-
eters discussed in previous sections. These assets range from fixed income products to
stock indices, and we summarize their return, risk, and correlations matrix in Figure 1.
Then we solve robust portfolio optimization problems using GI model and Lu model,
respectively. Since the only difference between Lu model and GI is in the uncertainty
set, we can make a meaningful comparison. Notice there is a hyper-parameter ω in
both models that represents the level of confidence, and it must be specified by user.
Consequently, we could have different optimal portfolio allocations based on different
values of ω.
 10

(a) Annualized return and risk of assets

(b) Correlation matrix of assets

Figure 1: Summary statistics of asset returns

Based on GI model, Portfolio allocations corresponding to different ω can be found

in Figure 2(a). When ω = 0, each uncertainty set shrinks to a single point, and the
robust optimization problem produces the same result as Markowitz’s model. The op-
timal portfolio says one should invest in 5 out of 13 assets. As ω increases, GI model
suggest that one allocate more in the first and sixth asset, and allocate less in the sec-
ond and third asset. Eventually, we see a 78% concentration in the 1st asset. In terms
of diversification, GI model doesn’t produce superior result in out study. One possible
explanation is uncertainty sets in GI model is too large due the projection methods used
in the model. As a result, the model tend to prefer few conservative assets with better
risk and return trade-off.

Portfolio allocations based on Lu model is shown in Figure 2(b). Notice the portfolios
displayed correspond to different value of c(ω), which is the critical value for ω based on
the distribution of joint ellipsoidal uncertainty sets. Since the distribution is difficult to
sample and optimized portfolios are not sensitive to changing ω, we use c(ω) as a proxy
of ω(see Section 2.5 for more on this issue). As c(ω) increases, we see similar trend as
in GI model, more allocation in the first and sixth asset and less in the second and the
third. However, Lu model picks up the seventh and the eleventh asset before completely
 11
abandoning the third asset and put over 70% weight in the first asset. Based on our
observation, Lu model does give slightly better diversification than GI model. The
observation is consistent with the claim that Lu model gives more accurate description
of the uncertainty sets.

(a) Portfolio allocations for different ω using GI model

(b) Portfolio allocations for different c(ω) using Lu model

Figure 2: Portfolio allocations using GI and Lu models

2.5 Discussion
To deal with statistical error and parameter perturbation, GI and Lu optimized the
objective function over uncertainty sets constructed around the least-square estimation
of µ and V. Because the objective is to optimize outcomes in worst scenarios, the
optimal solutions tend to be conservative. We observed quite conservative portfolio
allocations in numerical studies using both real market data and randomly generated
data. In addition to the conservativeness, numerical issues came up in our experimental
studies of these methods. For example, the inverse CDF of random variable Y = ni=1 Yi
P

in Lu model, which is the sum of n i.i.d. F −distributed random variable, is extremely

steep near the origin. As a result, the uncertainty set Sm,v (ω) and the optimal solution
is not sensitive to the confidence threshold ω, which is counterproductive to the primary
objective of robust optimization. To put these robust optimization models in practice,
 12
we need a more accurate description of underlying uncertainty sets and better numerical
technique to simulate complicated distributions.
Chapter 3

Bayesian Methods

The classic approach to portfolio optimization relies on estimations of market param-

eters such as mean return and covariance matrix. Yet many assumptions about the
distribution of asset returns are not realistic. Moreover, statistical methods used to es-
timate these parameters can be problematic. For example, estimating covariance matrix
when the sample size is small compared to the number of assets. As the optimization
problem is sensitive to input parameters, the resulting portfolio allocation can be sub-
optimal[24].

Bayesian theory provides a way to reduce the sensitivity of optimal portfolios to

input parameters. In Bayesian framework, we are able to consider not only statistical
estimations based on historical data, but also to incorporate prior information such as
investors’ subjective option or expert knowledge. By adjusting statistical estimations
with investors’ views, more reasonable and well-diversified portfolio allocations can be
achieved. In following sections, we focus on two models. The Black-Litterman model
was prominent for applying Bayesian method to estimation mean return and covariance
matrix. The robust Beyesian allocation approach extends the Black-Litterman model
by adding robust optimization elements to it.

13
 14
3.1 The Black-Litterman Model
Black and Litterman (BL)[7, 6] introduced a portfolio optimization framework which
incorporates more inputs than Markowitz’s model. Besides statistical estimation based
on market data, an investor is allowed to use her views on asset returns as input to the
underlying optimization problem. Consequently, portfolio allocations produced by BL
model naturally reflect investors’ opinion[32] in a consistent manner and do not need
post-processing afterwards.

Consider a market of N ∈ N securities or asset classes, whose returns X ∈ RN are

normally distributed with mean µ and covariance matrix Σ, i.e. X ∼ N (µ, Σ). Similar
to Markowitz’s model, Σ is estimated from exponential smoothing of historical time
series. However, BL treat µ as a random variable since its exact value is unknown. In
addition, BL assume distribution of µ is also normal with µ ∼ N (π, τ Σ) where τ ∈ R
is a small constant. Intuitively, π stands for the best estimation of µ and τ Σ represents
the uncertainty on this guess. Note the covariance matrix of µ is exactly the covariance
matrix of X scaled by the constant τ , reflecting the fact that µ is less volatile than X.

If there is no estimation error, then µ = π and X ∼ N (π, Σ). In the mean-variance

framework, an investor tries to find the optimal portfolio allocation from solving the
optimization problem
ω λ = argmax ω T π − λω T Σω
ω∈RN (19)
T
s.t. 1 ω = 1,
where the constant λ ≥ 0 represents the risk-aversion level which depends on each
individual investor. Assume there is an average risk-aversion level λ that applies to the
whole market, the quadratic optimization problem can be solved analytically. According
to the first order condition, the resulting market optimal portfolio ω
e satisfies

π = 2λΣe
ω (20)

and π is called the market equilibrium estimation. Thus we can set π in terms of
ω
e once the exogenous parameter λ is chosen (BL set it to be around 1.2). It is worth
pointing out that historical return information does not play a direct role in determining
 15
π. Assume historical returns {Xt }t=T
t=1 in different time periods are independent and
identically distributed (with the common distribution N (µ, Σ), the standard sample
estimator of µ is
t=T
1X 1
µ̂ = Xt ∼ N (µ, Σ), (21)
T T
t=1
1
Above equation suggests that it is reasonable to set τ = T. As discussed in previous
chapters, the sample estimator µ̂ tends to be inaccurate and its perturbation could lead
to dramatically different optimal portfolio in Markowitz’s model. One way to deal with
this issue is to apply Stein’s shrinkage estimator and set µ as a convex combination of
µ̂ and market equilibrium estimation π:

µs = (1 − s)µ̂ + sπ, (22)

where 0 ≤ s ≤ 1.

BL also allows investors to express their views on the expected return µ. An investor
can specify a “pick” matrix P ∈ RK×N , whose kth row describes a linear relative
performance relationship among assets. In addition, an investor can express level of
uncertainty in any of her views through a normal model

Pµ ∼ N (v, Ω). (23)

In short, the mean vector v ∈ RK describes relative performance among assets, and
the covariance matrix Ω quantifies levels of uncertainty of views. If an investor has a
specific view on asset returns, components of v can be chosen accordingly. However,
if an investor has only qualitative view, BL suggest v can be set in terms of market
volatility: q
vk = (Pπ)k + ηk (PΣPT )k,k , k = 1, ...., K, (24)

where ηk ∈ {−β, −α, α, β}. Here α and β are positive constants and different values of
ηk represent “bullish” or “bearish” views to some extent. Meucci[31] suggests that for
convenience one can set
1
Ω = PΣPT (25)
c
among other choices, where c is a positive constant representing investor’s confidence.
 16
Using Bayes formula, given P, v, and Ω, the posterior distribution of µ can be
obtained:
µ|v, Ω ∼ N (µBL , ΣBL ), (26)

where
µBL = π + τ ΣPT (τ PΣPT + Ω)−1 (v − Pπ)
(27)
ΣBL = (1 + τ )Σ − τ 2 ΣPT (τ PΣPT + Ω)−1 PΣ
d
The distribution of X can be computed using the fact that X = µ + Z where Z ∼
N (0, Σ). The posterior distribution of X is:

X|v, Ω ∼ N (µBL , ΣBL + Σ). (28)

Thus, an investor can use the posterior estimations µBL and (ΣBL + Σ) as the expected
return and covariance matrix in an mean-variance optimization model. BL compared
the resultant portfolios based on posterior estimations displayed meaningful differences
from the market portfolio, while preserving overall allocation characteristics.

3.2 Extension of Black-Litterman Model

Meucci considers how does an investor’s confidence in her views affect posterior parame-
ters µBL and ΣBL . He argues that when the confidence is full, the posterior distribution
N (µBL , ΣBL ) should converge the posterior distribution of N (π, Σ) conditioned on the
specific views v, that is
lim µBL = µ|v
Ω→0
(29)
lim ΣBL = Σ|v,
Ω→0

where

µ|v = π + ΣPT (PΣPT )−1 (v − Pπ)

(30)
Σ|v = Σ − ΣPT (PΣPT )−1 PΣ.
However, above relations are not completed satisfied in the original BL model. To
make the model more consistent with investors’ intuition, Meucci modifies BL model by
directly considering views on the asset return vector X ∈ RN . Instead of considering
 17
µ as a random variable, Meucci set it to be the market equilibrium estimation, i.e.
X ∼ N (π, Σ). In addition, the view is modelled as a normally distributed random
variable V with mean PX and covariance Ω. In particular, given any realization x of
X, the conditional distribution of V is

V|x ∼ N (Px, Ω). (31)

By Bayesian formula, the posterior distribution of X given a realization v of V and Ω

is
X|v, Ω ∼ N (µm m
BL , ΣBL ), (32)

where
−1
µm T T
BL = π + ΣP (PΣP + Ω) (v − Pπ)
(33)
−1
Σm T T
BL = Σ − ΣP (PΣP + Ω) PΣ.

Notice parameter τ in original BL model doesn’t appear in Meucci’s modification and

it’s no longer necessary to add the original covariance matrix Σ to the posterior. The
modified model can now perfectly match both the market and investors’ view. On the
m
one hand if Ω → ∞(no confidence in the view), the posterior distribution N (µm
BL , ΣBL )
converges to the prior distribution N (π, Σ). On the other hand if Ω → 0(full confi-
dence in the view), the posterior distribution converges to N (µ|v, Σ|v) that represents
investors’ view.

3.3 Robust Bayesian Allocation

So far, we have seen two approaches that aim to deal with statistical errors and port-
folios’ sensitivity. Robust optimization works by optimizing the worst-case outcome
over a set of possible values centered around estimations. Bayesian methods work by
adjusting statistical estimations toward investors’ views on asset returns. Meucci[31]
introduced a Robust Bayesian Allocation (RBA) model that combines robust optimiza-
tion and Bayesian methods. In RBA model, Bayesian methods are used not only to
combine the view and market data, but also to identify suitable uncertainty ranges for
the parameters such as mean return and covariance matrix. Then the robust optimiza-
tion problem is solved via convex optimization.
 18

Consider a market with N assets where trading happens at a set of discrete times.
An equivalent version of the classical mean-variance portfolio optimization problem is

max φT µ
φ∈RN

s.t. φT Σφ ≤ β, (34)

1T φ = 1,

where φ ∈ RN , µ ∈ RN , Σ ∈ RN ×N , and β are portfolio weights, mean return vector,

covariance matrix, and maximal acceptable portfolio variance. Note this problem is the
dual of Markowitz’s formulation defined in chapter 1. A robust version of the problem
can be set up similarly by considering sets of possible values for µ and Σ:

max min φT µ
φ∈RN µ∈Θµ

s.t. max φT Σφ ≤ β, (35)

Σ∈ΘΣ

1T φ = 1,

where Θµ and ΘΣ denote corresponding uncertainty sets for µ and Σ.

Assume at time t ∈ N an asset purchased τ ∈ N period ago is being traded at price

Pt
Pt ∈ R. We can define the return of this asset between time t−τ and t as Rt,τ = Pt−τ −1.
In addition, denote the return data of all assets between time t − τ and t as Rt,τ ∈ Rn .
In many linear factor models, Rt,τ is assumed to be normally distributed with mean µ
and covariance matrix Σ, i.e. Rt,τ |µ, Σ ∼ N (µ, Σ). Meucci models that the investor’s
prior knowledge about µ and Σ with a normal-inverse-Wishart distribution:
Σ
µ|Σ ∼ N (µ0 , )
T0
(36)
−1 Σ−1
Σ ∼ W(ν0 , 0 ),
ν0
where µ0 and Σ0 represent the investor’s view about µ and Σ. T0 and ν0 represent
the confidence in the view (in number of periods). The larger T0 ( or ν0 ) is the more
confidence an investor has in µ0 (or Σ0 ). The posterior knowledge on (µ, Σ) can be
updated with sample mean µ̂ and sample covariance Σ̂, which are estimated based on
 19
T ∈ N periods of historical market data. In Bayesian framework, there are analytic
formulas for the posterior estimation (µ1 , Σ1 ) of (µ, Σ):

T1 = T0 + T
ν1 = ν0 + T
T0 µ0 + T µ̂ (37)
µ1 =
T1
1 (µ − µ̂)(µ0 − µ̂)T
Σ1 = [ν0 Σ0 + T Σ̂ + 0 1 1 ].
ν1 T + T 0

In above equation, T1 and ν1 can be understood as the total confidence in terms of

the number of periods, and the posterior estimation (µ1 , Σ1 ) is the weighted average of
prior(µ0 , Σ0 ) and evidence (µ̂, Σ̂).

Meucci’s formulation not only gives the posterior estimation of (µ, Σ) but also the
uncertainty sets Θµ and ΘΣ :

Θµ = {µ : (µ − µce )T S−1 2
µ (µ − µce ) ≤ qµ }
(38)
ΘΣ = {Σ : vech[Σ − Σce ]T S−1 2
Σ vech[Σ − Σce ] ≤ qΣ },

where
µ̂ce = µ1
ν1
Σ̂ce = Σ1
ν1 + N + 1
1 ν1 (39)
Sµ = Σ1
T1 ν1 − 2
2ν12
SΣ = (DT (Σ−1 ⊗ Σ−1 −1
1 )DN ) ,
(ν1 + N + 1)3 N 1
where “vech” means stack the columns of a matrix skipping the redundant entries above
the diagonal, DN is the duplication matrix that reinstates the redundant entries above
the diagonal of a symmetric matrix, ⊗ is the Kronecker product, and (qµ , qΣ ) is a pair of
hyper parameters which determine the size of uncertainty sets. In short, the uncertainty
sets are ellipsoids centered at posterior estimations µ1 and Σ1 , which is quite similar
to their counterparts in the GI model.
 20
Meucci further proves that the robust optimization problem in equation (35) with
uncertainty sets define in equations (37) and (38) is equivalent to
q
T
max φ µ1 − γµ φT Σ1 φ
φ∈RN
(40)
s.t. φT Σ1 φ ≤ γΣ
1T φ = 1,

where r
qµ ν1
γµ =
T1 ν1 − 2
β (41)
γΣ = r .
ν 2ν12 qΣ
2

ν+N +1 + (ν+N +1)3

Finally, the optimization problem in equation (40) can treated as a SOCP and solved
numerically using existing convex optimization softwares such as CVX[17].

3.4 Numerical Study

The RBA model combines features of Bayesian methods and robust optimization. Con-
sequently, it requires input of both types of models we have previously discussed. For
example, prior mean return µ0 and prior covariance matrix Σ0 , together with their con-
fidence T0 and ν0 expressed in number of periods. In addition, we must also choose qµ
and qΣ that govern the size of uncertainty set (similar to the parameter ω is GI model).
Last but not least, another pair of mean and covariance-(µ̂, Σ̂)-needs to be estimated
from historical price data so that the prior estimation can be updated. We present some
summary statistics in Figure 3(a) including prior belief, estimations from data, and pos-
terior adjusted values of monthly return and risk of assets, which are the same as in the
study of GI or Lu models. In our numerical study, we assign equal confidence (weight)
to parameters that represent investor’s view and those estimated from historical data.
By using different levels of qµ and qΣ that control the size of uncertainty sets, we study
how portfolio allocations change accordingly.
 21

(a) Return and risk representing prior, data, and posterior estimations

(b) Portfolio allocations using RBA model

Figure 3: Summary statistics and portfolio allocation using RBA model

Portfolio allocations corresponding to different robustness (equal to qµ and qΣ ) can

be found in Figure 3(b). Overall, optimal portfolios are concentrated in 5 out of 13
assets. This might be caused by the risk-return feature of assets, and the fact that there
is no allocation constraint. In particular, when the robustness is zero, there is more
50% allocation in the ninth asset. When the robustness increases, the optimal portfolio
gradually picks up allocation in the fifth asset, while giving less weight on the other 4
assets. Eventually, the fifth asset get a 30% weight and the weight on the ninth asset is
reduced to be less than 50%. Within this limited case, RBA model does gives slightly
better diversification comparing to the classic model. However, we must emphasize
how much more complicated RBA is than Markowitz’s, or even the Black-Litterman’s
model. With many parameters that must be chosen by the investor, RBA definitely
requires more than average investment expertise. Together with sophisticated elements
in robust optimization, the model complexity is the greatest among all models we have
discussed so far. In summary, one should be very careful and have sufficient confidence
to use RBA model in practice.
Chapter 4

Dynamic Portfolio Optimization

So far, we have considered portfolio optimization in the static setting, where an investor
makes investment decision at a specific moment in time without further update. In
practice, new data and information are coming in every day, hour, and second. Statisti-
cal estimations and investors’ views can change frequently. Thus, we can treat portfolio
optimization as dynamic decision making problems. In this chapter we discuss two
dynamic programming problems. The first is a general robust dynamic programming
problem in discrete time. The second is a portfolio optimization problem which was
first studied by Merton.

In the sequential setting, investment decisions are made in different stages in time.
Each decision, besides generating a possible immediate reward, changes the ground of
all future decisions and rewards associated. Due to the uncertain nature of financial
markets, an investor only has limited information about the immediate reward of each
decision and future decisions are dependent on previous ones. To achieve a good per-
formance through all stages, the investor has to consider both the payoff of immediate
decision and how is this decision going to influence future payoffs. Dynamic program-
ming (DP) is a mathematical framework that allows the decision maker to compute
a good overall strategy by succinctly encoding the evolving information state. It was
originally developed by Richard Bellman[3] in the 1950s and has been widely applied in
areas such as aerospace engineering, computer science, and economics. The fundamental
idea behind DP is to break a complicated problem into nested sub-problems and solve

22
 23
them recursively. In formal DP theory, uncertainty in the environment is modeled by a
Markov process whose transition probability depends both on the information state and
the action taken by the decision maker. It assumes the transition probability of each
state-action pair is known exactly, and the goal is to choose a policy to maximizes an
objective function. The idea of dynamic programming can be extended to continuous
time with moderate mathematical sophistication. As we will see in Merton’s portfolio
problem, using stochastic calculus and differential equation, a continuous-time dynamic
programming can be elegantly formulated as a partial differential equation.

4.1 Robust Dynamic Programming

Iyengar[19] proposed a Robust Dynamic Programming (RDP) framework, where the
transition probability of each state-action pair is no longer assumed to be known ex-
actly (compare to classic dynamic programming problems).

Assume decisions are made at discrete point in time t ∈ T = {0, 1, ..., N − 1} where
N is a fixed positive integer. At each t ∈ T , the decision maker observes a current state
st ∈ St , where St is a discrete set of possible states at time t. Based the state st at time
t, the decision maker needs to take an action a ∈ A(st ) - also a discrete set of all possible
actions at st . For any discrete set B, let M (B) be the set of all probability measures
on B. Then a random action based on state st at time t naturally corresponds to a
probability measure qst ∈ M (A(st )). For t < N , the decision maker receives a reward
rt (st , at , st+1 ) if the action at is taken in state st at time t, and the state st+1 is reached
at time t + 1. The reward at the last stage rN is assumed to be a function of of sN only.

To describe the concept that the transient probability is not known exactly, Iyengar
introduced the notion of ambiguity. At each time t ∈ T , a set of probability measures
Pt (s, a) ⊂ M (St+1 ) is assigned to a state-action pair (s, a). This means at time t, if the
decision maker in state s takes action a, then the transient probability to all possible
states at t + 1 is an element of psa ∈ Pt (s, a). However, there is no way to determine
which element of Pt (s, a) is actually governing the transition process.
 24

Under the ambiguity of transition probability, a decision maker needs a rule to choose
suitable actions based on historical states and actions, denoted as ht = (s0 , a0 , ..., st−1 ,
at−1 , st ) for t < N . Let Ht be the set of all possible historical states and actions up
to time t. A randomized decision rule is a map dt : Ht 7→ M (A(st )), from Ht to the
set of probability measures on possible action based on state st . A policy π of a robust
dynamic programming problem is a collection of decision rule, i.e. π = {dt |t ∈ T }.
Given a decision rule dt , the set of probability measures consistent with dt is defined as

F dt = {P : Ht 7→ M (A(st ) × St+1 )|∀h ∈ Ht , Ph (a, s) = qdt (h) (a)pst a (s),

(42)
pst a ∈ P (st , a), a ∈ A(st ), s ∈ St+1 }.

It is also assumed (for regularity) that the set of probability measures consistent with a
policy π is just the product of sets of probability measures consistent with decision rule
dt at each time t < N :
F π = F d0 × F d1 × ... × F dN −1 . (43)

In a RDP problem, the reward V0π (s) generated by a policy π starting from the
initial state s0 = s is defined as
X
V0π (s) = inf π EP [rN (sN ) + rt (st , dt (ht ), st+1 )], (44)
P ∈F
0≤t<N

where EP denotes the expectation w.r.t. the probability measure P ∈ F π . The in-
terpretation is that the reward of a policy π is the minimal expected reward over all
probability measures consistent with π. The goal of a decision maker is to find an
optimal policy π ∗ that maximize the robust value function

V0∗ (s) = sup V0π (s), (45)

π∈Π

where Π is the set of all admissible history-dependent policies. Similarly, the reward of
policy π starting from history hn with n < N is defined as
X
Vnπ (hn ) = inf π EP [rN (sN ) + rt (st , dt (ht ), st+1 )], (46)
P ∈Fn
n≤t<N
 25
where Fnπ = F dn × F dn+1 × ... × F dN −1 . The optimal reward starting from hn is defined
as
Vn∗ (hn ) = sup Vnπ (hn ), (47)
π∈Πn

where Πn denotes the set of all admissible history-dependent policies concerning t ≥ n.

One remarkable result on RDP is the robust Bellman equation. It states that optimal
reward functions {Vn∗ |0 ≤ n ≤ N } satisfy the recursive relation

VN∗ (hN ) = rN (sN ),

(48)
Vn∗ (hn ) = sup { inf ∗
Ep [rn (sn , a, s) + Vn+1 (hn , a, s)]}, n < N,
a∈A(sn ) p∈P (sn ,a)

∗ (h , a, s) means the optimal reward function at time n + 1 with history date

where Vn+1 n

obtained by concatenating a and s to hn . The relation suggests that the optimal reward
given history hn can be achieved by considering all possible action in A(sn ) and move
one step forward in time. Although the relation gives a guideline of the connection
between optimal reward functions at different times, it is not easy to apply due to fact
that every state and action in all previous times come in to play. Fortunately, there is
a significantly simplified equivalent version which states the following:

VN∗ (sN ) = rN (sN ),

(49)
Vn∗ (sn ) = sup { inf ∗
Ep [rn (sn , a, s) + Vn+1 (s)]}, n < N,
a∈A(sn ) p∈P (sn ,a)

where s ∈ Sn+1 . Compare to the history-dependent version, above relation says robust
value functions Vn∗ (hn ) are functions of the most recent states only, i.e. in the form of
Vn∗ (sn ). Instead of carrying a sequence of historical events, we only need to consider
the current state and the immediate future. Thus, the robust dynamic programming is
a Markov decision-making problem with ambiguity in transition probability.

When dealing with discrete states and actions, RDP can be a powerful tool. How-
ever, it is not easy to find a good, natural set of ambiguity sets. In addition, it’s
not completely clear how to adapt RDP to continuous decision-making problems. Al-
though approximation by discretizing time, action, and underlying probability measures
 26
is a possibility, the computational complexity and numerical issues that come with dis-
cretization can be difficult to deal with.

4.2 Merton’s portfolio problem

Robert Merton[28, 30] examined the problem of optimal portfolio selection and con-
sumption rules in continuous-time setting. Unlike discrete-time portfolio allocation
problems, assets continuously generate returns in Merton’s model and they are the only
source of an investor’s income. Likewise, consumption and portfolio allocations can also
be continuously adjusted by the investor.

Assume there are m ∈ N assets in the market and the price of ith asset, Pi (t)
is generated generated by a Geometric Brownian Motion (GBM) with constant drift
ai ∈ R and constant volatility σi ∈ R, i.e.

dPi (t) = ai Pi (t)dt + σi Pi (t)dBi (t), (50)

where {Bi (t)}1≤i≤m are independent Brownian Motions (BM). At time t ∈ R, denote the
amount of total wealth as W (t), consumption per unit time as C(t), and the proportion
of total wealth in the ith asset as ωi (t) for 1 ≤ i ≤ m. It is understood that the investor
distribute all wealth among m assets, i.e. m
P
i=1 ωi (t) = 1 for any t. It is easy to derive a
Stochastic Differential Equation (SDE) that describes the evolution of W (t) over time:
"m # m
X X
dW (t) = ωi (t)ai W (t) − C(t) dt + ωi (t)σi W (t)dBi (t). (51)
i=1 i=1

Merton further studied a special case of above problem where there are only two
assets in the market. One is a risky asset with drift a ∈ R and volatility σ > 0. The
other is risk-free with drift r ∈ R and zero volatility. In the two-asset scenario, if we
denote the proportion of total wealth in the risky asset as ω(t), then the proportion in
the risk-free asset is 1 − ω(t). The governing equation for the wealth process becomes

dW (t) = [(r + (a − r)ω(t))W (t) − C(t)]dt + ω(t)σW (t)dB(t), (52)

 27
where B(t) is a standard Brownian Motion (BM). Given the equation, an investor’s
goal is to find the optimal consumption rule C(t) and portfolio allocation ω(t) which
maximize the sum of expected utility and bequest at some future time:
Z T
max E[Q(W (T ), T ) + e−ρτ U (C(τ ))dτ ]
C(t),ω(t) 0

s.t. C(t) ≥ 0 (53)

ω(t) ≥ 0
W (0) = W0 ≥ 0.

In above formulation, E is short for the the conditional expectation operator E(0) at ini-
tial time 0, T stands for the the investment horizon, ρ is a non-negative discount factor,
Q(·, ·) is the bequest valuation function, and U (·) is the utility function. The first part
inside the expectation represent the bequest at time T and the second part stands for the
accumulated utility incurred by consumption between time 0 and T . In addition, stan-
dard utility theory assumes U (·) is concave and smooth enough with U 0 > 0 and U 00 < 0.

To derive the optimal conditions for the weight ω(t) and consumption C(t), Merton
utilized the Bellman principle of optimality for Dynamic Programming. For t ∈ [0, T ],
define Z T
I(W (t), t) = max E(t)[ e−ρτ U (C(τ ))dτ + Q(W (T ), T )]
C(s),ω(s) t
(54)
s.t. C(s) ≥ 0
ω(s) ≥ 0,
where E(t) is the conditional expectation operator at time t. I(W (t), t) is understood
as the maximal expected utility that can be achieved under any consumption C(s)
and weight ω(s) starting at time t. The principle of DP implies that for any t0 < t,
I(W (t0 ), t0 ) should be determined by the following optimization problem
Z t
I(W (t0 ), t0 ) = max E(t0 )[ e−ρτ U (C(τ ))dτ + I(W (t), t)]
C(s),ω(s) t0
(55)
s.t. C(s) ≥ 0
ω(s) ≥ 0.
 28
Assume I(·, ·) has bounded third partial derivatives, we can apply Taylor’s theorem to
I(W (t), t) with t = t0 +h and let h →
− 0, the governing equation for optimal consumption
C(t) and ω(t) at any t ∈ [0, T ] can be derived:

∂I ∂I
0 = max {e−ρt U (C(t)) + + (ω(t)[(a − r) + r)W (t) − C(t)]
C(t),ω(t) ∂t ∂W
(56)
1 ∂2I 2
+ ω (t)σ 2 W 2 (t)}.
2 ∂W 2
Denote the objective function to be maximized in the above equation as Φ(ω, C, W, t).
At a fixed time t, the first-order conditions for a regular interior maximum are

 ∂Φ = 0 = e−ρt U 0 (C) − ∂I
∂C ∂W
(57)
∂Φ ∂I ∂2I

∂ω =0= ∂W (a − r) + ∂W 2
ωσ 2 W.

∂2Φ ∂2Φ ∂2Φ

For the second order conditions, note ∂C∂ω = ∂ω∂C = 0, and ∂C 2
= e−ρt U 00 (C) < 0
since U is assumed to be strictly concave. Merton further assumed I(W, t) is strictly
concave in W . Consequently, the second order conditions are automatically satisfied.
Thus, the original problem of determining optimal consumption and portfolio allocation
is reduced to solving the following system of equations



 Φ(ω ∗ , C ∗ , W, t) = 0

Φω (ω ∗ , C ∗ , W, t) = 0 (58)


ΦC (ω ∗ , C ∗ , W, t) = 0



under the boundary condition I(W (T ), T ) = Q(W (T ), T ) at the terminal time.

In general, the system of equations in (58) are difficult to solve due to the coupling
between partial differential equation and algebraic equations. However, it can be solved
explicitly if the utility has constant relative risk-aversion, i.e.
Cγ
U (C) = , (59)
γ

where γ ∈ (0, 1). If we further assume that the bequest at the terminal time is given by
 29
Q(W (T ), T ) = 1−γ e−ρT U (W (T ))(0 <   1), the optimality conditions become
 1


 C ∗ (t) = [eρt ∂W
∂It 1−γ
]
 ∂I
−(a−r) t

ω ∗ (t) = 2 2 ∂ ∂W

σ W
2I
t (60)
∂W 2

−ρt h i γ ∂I
(a−r)2 [ ∂Wt ]2

1−γ ∂It γ−1
+ ∂I ∂It

 0=e ∂t + ∂W rW − 2σ 2 ,

 1−γ t
γ ∂W ∂ 2 It
∂W 2

γ
subject to I(W (T ), T ) = 1−γ e−ρT [W (Tγ
)]
. To solve the last equation in (52), sub-
¯ b(t) −ρt
stitute a trial solution It (W (t, t)) = e [W (t)]γ and it turns out b(t) must satisfy
γ
the ordinary differential equation
−γ
b0 (t) = µb(t) − (1 − γ)[b(t)] 1−γ
(61)
s.t. b(T ) = 1−γ ,
h 2
i
where the constant µ = ρ−γ 2σ(a−r)2 (1−γ) + r . The governing equation of b(t) has solution
vt−T 1−γ
h i
b(t) = 1+(v−1)e
v , where v = 1−γ µ
. Based on b(t), the optimal solutions for ω(t)
and C(t) are
a−r
ω ∗ (t) =
σ 2 (1
− γ)

 v
W (t) if v 6= 0 (62)
1+(v−1)evt−T
C ∗ (t) =
1
T −t+ W (t) if v = 0.


There is a very important implication from the expression of ω ∗ (t). Notice ω ∗ (t) is a
constant independent of time. The solution to Merton’s problem suggest that with any
constant relative risk-aversion utility function, an investor should put a constant portion
of wealth into the risky asset to maximize the current expectation of utility at terminal
time. Notably, Merton also studied a wide range of related problems which involved
other types of utility function, multiple asset classes, and finite/infinite time horizon.
His pioneering work had tremendous impacts on portfolio allocation, asset pricing, and
the study of financial derivatives. In 1997, Robert Merton and Myron Scholes were
awarded the Nobel Memorial Prize in Economic Sciences for their work in derivative
pricing in collaboration with the late Fischer Black (1938-1995).

Nowadays, Merton problem has become a classic example in stochastic control[33]

and we do not need to go through all the technical steps went through by Merton. In
 30
general, there are mainly three steps to solve a stochastic control problem. The first
step is to transform the problem into a Partial Differential Equation (PDE) using Dy-
namic Programming Principle. The second step is to solve the PDE, usually by further
converting it to a system of Ordinary Differential Equations. Finally, if there is a nice
solution to the PDE, we come back and try to verify it is actually a solution to the
original stochastic control problem.
Chapter 5

Dynamic Portfolio Optimization

with Bayesian Update - Discrete
Model

We consider a discrete portfolio optimization problem that is in the intersection of

dynamic programming, robust optimization, and Bayesian theory. To be specific, an
investor can adjust portfolio allocation at the beginning of each period, which makes
the problem dynamic. However, a crucial parameter that governs the transitional prob-
ability measure of price movement is unknown to us. Consequently, we cannot solve the
problem directly using Dynamic Programming principle. To deal with the parameter
uncertainty, we model the unknown parameter via a probability distribution. Together
with Bayesian theory, we are able to update the distribution of unknown parameter
with price movement we have observed. Eventually, we will see how to solve it with
dynamic programming and Bayesian update.

5.1 A Discrete Model with an Uncertain Parameter

We consider a discrete-time market where there are T ∈ N investment periods (indexed
by 0, 1, ..., T −1). In this market, there is a risky asset and a risk-free asset. The risk-free

31
 32
asset offers a constant rate of return r ∈ R in each period. In contrast, the price of
risky asset can either go up by a factor ea or go down by a factor e−a , where a ∈ R+
is a constant. However, the probability of price increase is not known exactly. We
treat the probability of price increase as a continuous random variable and use the beta
distribution to model its probability distribution. The beta distribution[21] is a family
of continuous probability distributions defined on the interval [0, 1] and parameterized
by two positive shape parameters α ∈ R and β ∈ R. The probability density function
given α and β is
Γ(α + β) α−1
B(x, α, β) = x (1 − x)β−1 , (63)
Γ(α)Γ(β)
where Γ(·) is the gamma function[13]. Thus, the probability that the price of risky
asset will increase in a certain period could be any number in [0, 1] with respect to
the beta distribution B(x, α, β). We assign a Bayesian-type method for the investor to
make inference about the probability distribution. Initially, we set α = β = 1. For
0 ≤ t < T , denote it as the number of periods where the price of risky asset moves up
in the previous t periods, and jt as the number of periods where the price moves down.
A simple counting argument implies us it + jt = t. Based on it and jt , the investor then
update her belief about the distribution as B(x, 1 + it , 1 + jt ). Under this mechanism,
the investor believes the distribution is likely to tilts toward 1 if jt > it and vice versa.

With the price dynamics and belief update mechanism laid out, we consider the
investment problem. At the beginning of each period 0 ≤ t < T , the investor decides
to put πt ∈ [0, 1] portion of her wealth in the risky asset, and rest 1 − πt portion goes
to the risk-free asset. The investor then hold the portfolio until the beginning of next
period when a reallocation is allowed. Depending on the previous decision and actual
price movement, the value of portfolio could go up or down. Similar to the Merton
problem, the objective of the investor is to find the optimal portfolio allocation πt∗ so
that the expected utility of at some future time is maximized. In our problem, we do
not consider consumption of portfolio. The utility at the end of the T th period is a
strictly concave function of portfolio value. For any time period t < T , the utility is
determined by maximizing the expected utility of immediate future over all possible
investment decisions. Formally, define the utility function at the end of the T th period
 33
as
y 1−γ
UT (y) = , (64)
1−γ
where y denotes the total portfolio value at the end of T th period, and γ 6= 1 is a positive
constant. For 0 ≤ t < T , denote the portfolio value at the beginning of tth period as
vt ∈ R. If the investor decide to put πt ∈ [0, 1] portion of vt in the risky asset, her
portfolio value at the beginning of next period (or the end of this period, equivalently)
could be vt (1 − πt )er + vt πt ea if the price of risky asset increase during tth period, or
vt (1 − πt )er + vt πt e−a if the price decrease. Since the actual probability of price change
is only known as a distribution, we define the value function associated with portfolio
value vt and allocation πt as
Z 1 h
Vt (vt , πt ) = Ut+1 (vt (1 − πt )er + vt πt ea )x
0
i (65)
r −a
+ Ut+1 (vt (1 − πt )e + vt πt e )(1 − x) B(x, 1 + it , 1 + jt )dx.

In other words, Vt (vt , πt ) is the expected utility of immediate future w.r.t. the distribu-
tion of price move that represents investor’s current belief. Based on Vt (vt , πt ), we can
define the utility of vt as
Ut (vt ) = max Vt (vt , πt ), (66)
πt ∈[0,1]

which is the maximal value function among all possible allocations. The optimal port-
folio allocation is naturally given by

πt∗ = argmaxVt (vt , πt ). (67)

πt ∈[0,1]

In summary, the investor chooses πt∗ that maximize the value function of vt for each
0 ≤ t < T based on her current belief of the probability distribution. In particular, at
the very beginning of the investment process she wants to find

π0∗ = argmaxV0 (v0 , π0 ). (68)

π0 ∈[0,1]
 34
5.2 Transformation to a Dynamic Programming
It turns out the portfolio allocation problem described in section 5.1 can be reformulated
into an ordinary dynamic programming problem. We describe details in the following
theorems.

Theorem 1. For 0 ≤ t ≤ T , utility function Ut (vt ) takes the form

Ut (vt ) = ct (it , jt )vt1−γ , (69)

where the coefficient ct is a function of it and jt which govern the beta distribution
B(x, 1 + it , 1 + jt ).

Proof. We proceed by backward induction. First note the utility function after the last
y 1−γ
investment period has the form UT (y) = 1−γ , which satisfy the claim. Now assume
1−γ
Un+1 (vn+1 ) = cn+1 (in+1 , jn+1 )vn+1 satisfy the claim for some n ≤ T − 1. Consider the
value function for any pair (vn , πn )
Z 1h
Vn (vn , πn ) = Un+1 (vn (1 − πn )er + vn πn ea )x
0
i (70)
r −a
+ Un+1 (vn (1 − πn )e + vn πn e )(1 − x) B(x, 1 + in , 1 + jn )dx.
By the inductive assumption about Un+1 , we can rewrite the above integral as
Z 1h
Vn (vn , πn ) = cn+1 (in + 1, jn )vn1−γ (er − πn er + πn ea )1−γ x
0
i
+ cn+1 (in , jn + 1)vn1−γ (er − πn er + πn e−a )1−γ (1 − x) B(x, 1 + in , 1 + jn )dx,
(71)
where the difference in cn+1 reflects different price movement in two different scenarios.
After pulling vn1−γ outside, the integral is then completely determined by in , jn , and
πn . Note when evaluating the integral, in and jn have been determined from past price
movement. Thus we can write

Vn (vn , πn ) = fn (in , jn , πn )vn1−γ , (72)

where
Z 1h
fn (in , jn , πn ) = cn+1 (in + 1, jn )(er − πn er + πn ea )1−γ x
0
i
+ cn+1 (in , jn + 1)(er − πn er + πn e−a )1−γ (1 − x) B(x, 1 + in , 1 + jn )dx.
(73)
 35
Now the utility function can be written as

Un (vn ) = vn1−γ max fn (in , jn , πn ), (74)

πn ∈[0,1]

where the last term depend only on in and jn by maximizing over πn . Thus the inductive
step is complete and claim is valid for all 0 ≤ t ≤ T .

Consequently, the theorem suggest that the amount of portfolio value vt does not
affect the optimal allocation, which is consistent with our intuition since we talk about
portion πt rather than amount. Moreover, it is the historical price movement it and jt
that determines the optima allocation through the beta distribution distribution. After
separating vt , at each state (it , jt ), we only need to solve a one dimensional optimization
problem over πt ∈ [0, 1]. However, the expression of fn (in , jn , πn ) makes it very difficult
to find a closed-form solution to optimal, mainly due to the fractional exponents and
the fact we restrict πt within [0, 1]. Fortunately, in practice we don’t need to consider
every possible allocation between 0 and 1, and a discrete approximation should serve
the purpose, e.g. all integer multiples of 1% or 0.1%.

When πt is restricted to be chosen from a finite action set A, the maximization

problem can be solved numerically by enumerating all choices of πt . At time t, the
portfolio optimization problem with the pair (it , jt ) becomes

argmaxft (it , jt , πt ), (75)

πt ∈A

where ft (it , jt , πt ) is defined in the proof of the previous theorem. Combining the pro-
cess presented in the proof, we summarize the complete problem-solving process in the
algorithm below.
Algorithm 1: Discrete Dynamic Programming with Bayesian Update
Input: number of periods T , action set A,
1 For t = T − 1, T − 2, ...1, 0
(i) For every possible combination of it and jt , solve the discrete optimization
problem
argmaxft (it , jt , πt )
πt ∈A

to find the the optimal πt∗ and coefficient ct (it , jt ).

(ii) Stop if t = 0. Otherwise continue to next value of t.
 36
Note in above algorithm, we must solve the optimization problem for all possibilities
of it and jt for a given period, which requires ct+1 (it+1 , jt+1 ) at the next time period.
That is why we must solve the dynamic programming problem backward in time. Due
to the constraint it + jt = t, there are t + 1 different pairs of (it , jt ). The connections
between states (it , jt ) at different times can be conveniently represented by a binomial
tree, which can be efficiently modeled in many numerical computing software. We will
present some result of numerical experiment in the next section.

5.3 Numerical Study

We consider two instances of numerical studies: one does not allow short sell and the
other does. In both cases we set the total number of periods T = 8, the return of
risk-free asset r = 0.1, price jump constant a = 0.1, and the constant in utility func-
tion γ = 2. The dynamic programming algorithm described earlier is implemented in
Python.

In the case where no short sell is allowed, the admissible control set is
n n o
A= |n ∈ N, 0 ≤ n ≤ 1000 . (76)
1000
The result of numerical study is presented in Figure 4. In Figure 4, each branch rep-
resents a possible price movement and each node represents a point when the investor
must make decision about πt . The number attached to each branch is the expected
probability (according to the distribution B(x, it , jt )) of price moving up or down, de-
pending on whether the branch leads to a higher of lower node. If the branch is green
(or red), it leads to a higher (or lower) price in the next period. The probability of
price movement can also be visualized in the width of each branch. For example, if the
price keeps decreasing, the probability of having an even lower price in the next period
will increase, which is reflected by the increasing width of red branches as well as the
increasing values associated with them. The number attached to each node directs how
much (in percentage) the investor should put in the risky asset. The general pattern
is that if the price drops more than it moves up, the investor would avoid risking the
 37
portfolio at all. For example in period 7, the optimal allocation is to completely avoid
the risky asset if one were at 4 nodes from the bottom. On the other hand, if the price
moves up more often than it drops, the investor starts to invest in the risky asset, e.g.
4 nodes from the top in period 7. In general, this model leads to very conservative in-
vestment decision and the reason is that there is no “drift” term that drives the return
of risky assets, i.e. the expected return of risky asset has no obvious advantage over the
risk-free asset.

Figure 4: Optimal action tree (no short sell allowed)

Next we consider the situation where short sell is allowed. While everything else is
unchanged, the admissible control set becomes
n n o
A= |n ∈ Z, −1000 ≤ n ≤ 1000 , (77)
1000
which means an investor can choose to short the risky asset up to 100%. The result of
numerical study is presented in Figure 5. The general pattern is similar to what we can
 38
observe from Figure 5, i.e. optimal allocations depends how does price move in previous
periods. However, since short sell is allowed, it becomes beneficial to short when price
of the risky asset becomes really low, e.g. bottom nodes at each time period.

Figure 5: Optimal action tree (short sell allowed)

Chapter 6

Dynamic Portfolio Optimization

with Bayesian Update -
Continuous Model

Models we have discussed in previous chapters can be seen as different extensions

to Markowitz’s model. Merton’s portfolio problem studies portfolio optimization in
continuous-time where an investor can continuously adjust portfolio allocations. The ex-
pected return and covariance in Markowitz’s model correspond to the drift and volatility
of the underlying price process in Merton’s model. Bayesian methods allows statistical
estimations to be modified according to investor’s opinions or additional observations.
Consequently, model parameters can also be updated in a dynamic fashion, as we have
seen in discrete Bayesian models. Thus, it is natural to consider the combination of
Bayesian methods and continuous-time models. In this chapter, we shall develop a
continuous Bayesian model and study a portfolio optimization problem based on the
model.

39
 40
6.1 Conjugate Distributions
In Bayesian statistics, we often want to infer the distribution of a parameter θ based on
observed data x. According to Bayes’ formula[8]
p(x|θ)p(θ)
p(θ|x) = , (78)
p(x)
the posterior distribution p(θ|x) is determined by the product of prior distribution p(θ)
and likelihood function p(x|θ). If the posterior distributions p(θ|x) are in the same
probability distribution family as the prior probability distribution p(θ), the prior and
posterior are called conjugate distributions, and the prior is called a conjugate prior for
the likelihood function.

Let X be a normal random variable with distribution X ∼ N (µ, σ 2 ), where the variance
σ 2 is a known constant but the mean µ is not known for sure. Assume µ is also a normal
random variable with (prior) distribution µ ∼ N (µ0 , σ0 2 ). In the Bayesian framework,
we can update the (posterior) distribution of µ from independent observations of X.
If there are n ∈ N independent observations D = (x1 , ..., xn ) of X, then the likelihood
function is
n Pn
n 1 2
p(xi |µ) = (2πσ 2 )− 2 e− 2σ2 i=1 (xi −µ)
Y
p(D|µ) = . (79)
i=1
According to Bayes’s formula, the posterior probability density of µ is a constant mul-
tiple of Pn
1 2 1 2
− 2 (µ−µ0 ) − 2σ 2 i=1 (xi −µ)
2σ0
e . (80)

Notice the posterior distribution of µ is also normal with µ|D ∼ N (µp , σp2 ), where
Pn
xi
1
σp2
= 1
σ02
+ n
σ2
, and µp = σp2 ( σµ02 + i=1
σ2
). As we can see, the prior and posterior
0
distribution of µ are both Gaussian, hence conjugated distributions. With the posterior
distribution of µ computed, we can further update the posterior predictive distribution
of X by Z
p(x|D) = p(x|µ)p(µ|D)dµ
Z 1 2
(81)
1 1 2 1 − 2 (x−µp )
= (2πσ 2 )− 2 e− 2σ2 (x−µ) (2πσp2 )− 2 e 2σp
dµ.

A simple computation shows the posterior predictive distribution of X is again normal,

i.e. X|D ∼ N (µp , σp2 + σ 2 ).
 41
6.2 Stochastic Market Model
In this section, we will extend the idea of conjugate distributions to a continuous stochas-
tic process with unknown drift. Consider a stochastic process Xt evolves according to
the stochastic differential equation (SDE)

dXt = µdt + σdWt , (82)

where Wt is a standard one dimensional Brownian motion, σ is a positive constant, and

µ is a constant but its value is not known for sure. Based on information available at
time t = 0, our belief about µ is a normal distribution with mean µ0 and standard de-
viation σ0 , i.e. µ ∼ N (µ0 , σ02 ). This distribution represents our belief about µ without
observing Xt . Later when we observe values of Xt , the distribution of µ can be updated
using the conjugate technique described earlier.

For example, if we have observed Xt on the time interval [0, s], we can apply the
following approach to update the distribution of µ. First for a positive integer N , par-
s
tition [0, s] into N sub-interval with {si = i∆s|0 ≤ i ≤ N, ∆s = N }. Then define
random variables ∆Xi = Xsi − Xsi−1 for 1 ≤ i ≤ N . Due to the stationary indepen-
dent increment property of Brownian motion, we know {∆Xi }1≤i≤N are independent
and identically distributed with a common normal distribution N (µ∆s, σ 2 ∆s). By our
prior knowledge of µ, the random variable µ∆s also has a (prior) normal distribution
N (µ0 ∆s, σ02 ∆s2 ). Now we can repeat the conjugate distributions procedure to update
the posterior distribution of µ∆s(hence µ via scaling) given {∆Xi }1≤i≤N , that is
∆sµ0 σ 2 + ∆sσ02 Xs σ 2 σ02 (∆s)2
 
µ∆s|{∆Xi }1≤i≤N ∼ N , 2 , (83)
(σ 2 + sσ02 ) σ + sσ02
where we have used the identity N
P
i=1 ∆Xi = Xs . Similarly, the posterior predictive
distribution of the increment ∆X = Xs+∆s − Xs , becomes
∆sµ0 σ 2 + ∆sσ02 Xs σ 2 σ02 (∆s)2
 
2
∆X|{∆Xi }1≤i≤N ∼ N , 2 + σ ∆s . (84)
(σ 2 + sσ02 ) σ + sσ02
Finally let N →
− ∞ and normalize the posterior predictive distribution over ∆s, we
obtain a normal distribution with following mean and variance
∆sµ0 σ 2 +∆sσ02 Xt
(σ 2 +sσ02 ) µ0 σ 2 + σ02 Xs
lim = , (85)
∆s→−0 ∆s (σ 2 + sσ02 )
 42
and
σ 2 σ02 (∆s)2
σ 2 +sσ02
+ σ 2 ∆s
lim = σ2. (86)
∆s→
−0 ∆s
Above equations imply that if we update the distribution of µ using observations up to
time s, the posterior predictive distribution of the normalized increment of Xt starting
Xs+∆s −Xs µ0 σ 2 +σ02 Xs
at time s, i.e. ∆s , is N ( (σ 2 +sσ02 )
, σ 2 ). If we further draw an analogy between
the discrete updating process and the incremental distribution of the stochastic process
described by equation (82), we obtain the following SDE corresponding to the posterior
predictive distribution
µ0 σ 2 + σ02 Xt
dXt = dt + σdWt . (87)
(σ 2 + tσ02 )
Notice equation (82) cannot be used directly to set a stochastic control problem since the
value of µ is unknown. However, we were able to derive equation (87) in the Bayesian
framework which avoids direct use of µ by combining the prior distribution of µ and
observations of the stochastic process. Thus, it becomes possible to study stochastic
control problems based on equation (87). It is easy to verify
Z t
2 2 σ
Xt = µ0 t + (σ + tσ0 ) 2 2 dWs (88)
0 σ + sσ0

assuming X0 = 0, from which we see Xt is a Gaussian process with normal distribution

 2

N µ0 t, σσ2 (σ 2 + tσ02 )2 ( σ12 − σ2 +tσ
1
2 ) at any time t.
0 0

Remark. From equation (83), we can derive the posterior distribution of µ based on
observations up to time t, denoting as µt ,
µ0 σ 2 + σ02 Xt σ 2 σ02
 
µt ∼ N , . (89)
(σ 2 + tσ02 ) σ 2 + tσ02
In stochastic filtering theory, the conditional expectation of µ over the σ−algebra gener-
ated by observations of {Xs }0≤x≤t has the same distribution as in equation (89). For a
more rigorous derivation within the measure-theoretic framework, one can refer to Bain
and Crisan[2].

6.3 Stochastic Control Problem

Consider a market where there is one risky asset, and one risk-free asset. The price
process of the risky asset follows the stochastic process Yt = eXt , where Xt is defined in
 43
equation (87). By Itô’s lemma, we can derive a stochastic differential equation for Yt
1

dYt = f (t, Yt ) + σ 2 Yt dt + σYt dWt , (90)
2
where
µ0 σ 2 + σ02 log y
f (t, y) = . (91)
σ 2 + tσ02
Assume the risk-free asset provides a constant (continuously compounding)rate of return
r, which leads to the ordinary differential equation for its price process

dBt = rBt dt. (92)

Based on market dynamics, an investor can continuously decide the portion of her
money to put in the risky asset (the rest goes to the risk-free asset). Her objective is
to maximize the expected utility of pay-off at some future time T . Denote the total
value of her investment at time t as Pt , then Pt should satisfy the following stochastic
differential equation if αt portion of Pt were invested in the risky asset
dYt dBt
dPt = αt Pt + (1 − αt )Pt , (93)
Yt Bt
which is equivalent to
h 1 2 i
dPt = αt (f (t, Yt ) + σ − r) + r Pt dt + αt σPt dWt . (94)
2
We further set the investor’s pay-off at a fixed future time T to be a constant relative
risk aversion (CRRA) utility function of the portfolio value PT ∈ R+
PT1−γ
U (PT ) = , (95)
1−γ
where γ > 0 is a constant that is not equal to 1 (in this thesis we will mainly focus on
the case where γ > 1). For any t ≤ T the investor has the freedom to choose αt ∈ A
so that the expected pay-off can be maximized, where A is a pre-defined action set
which describes possible actions the investor can take. For example, A = R means
the investor can take arbitrary long/short positions in both assets, A = {0.5} means
the investment must split the portfolio evenly into two assets. Formally, define the
value function V (t, y, p) as the conditional expectation given the price of risky asset and
portfolio value at time t

V (t, y, p) = sup E[U (PT )|Yt = y, Pt = p], (96)

αt ∈A
 44
where 0 ≤ t ≤ T , y, p ∈ R+ .

The standard approach to solve a stochastic control problem can be decomposed into
several steps. The first step is to set up the Hamilton-Jacobi-Bellman (HJB) equation
that the value function should satisfy. The second step is to obtain a solution to the
HJB equation, e.g. by solving a system of ODE. The final step is to verify the solution to
the HJB equation is indeed the value function of the stochastic control problem. Under
some regularity conditions, the value function is a smooth solution to the HJB equaion,
for example Friedman[14], Gilbarg[15],and Krylov[23]. However, those conditions are
not met in our problem. Yet we will see a smooth solution can be found when T satisfies
some conditions.

The HJB equation of this stochastic control problem is

n 1 1
0 = Vt + sup Vy (f + σ 2 )y + Vp [αt (f + σ 2 − r) + r]p
αt ∈A 2 2
1 1 o
+ Vyy σ 2 y 2 + Vpp αt2 σ 2 p2 + Vyp αt σ 2 yp . (97)
2 2
The process of deriving HJB equation in a stochastic control problem typically involves
Itô’s Lemma, and we include details in the Appendix. In the coming sections, we will
look at three different cases where the action set A is the set of all real numbers, contains
a single real number, and is a closed interval. To make it notation-wise easier, we denote
the expression inside the curly bracket in equation (97) as

1 1
G(αt ) = Vy (f + σ 2 )y + Vp [αt (f + σ 2 − r) + r]p
2 2
1 1
+ Vyy σ 2 y 2 + Vpp αt2 σ 2 p2 + Vyp αt σ 2 yp, (98)
2 2
emphasizing that we need to solve a maximization problem within equation (97). Notice
that G(αt ) is a quadratic function of αt and we can obtain a closed-form expression in
some cases.
 45
6.4 Stochastic Control Problem - Unconstrained Action
We first look at the case where the action set A is the entire real line, i.e. A = R.
Clearly the function G attains its maximum at

Vp (r − f − 21 σ 2 ) − Vyp σ 2 y
αt∗ = . (99)
Vpp σ 2 p

Substitute αt∗ back to the HJB equation (97) yields a second order non-linear partial
differential equation
1 1 1 1
0 = Vt + Vy (f + σ 2 )y + prVp + Vyy σ 2 y 2 − [Vp (f + σ 2 − r) + Vyp yσ 2 ]2 . (100)
2 2 2Vpp σ 2 2

Notice the stochastic differential equation (94) is linear in the process Pt , meaning if
Pt solve the equation so does any constant multiple of Pt . Together with the utility
function (17) at time T , we can assume the value function V has the following form

V (t, y, p) = A(t, y)p1−γ . (101)

Substituting equation (101) into equation (100) gives us a non-linear partial differential
equation which does not involve p

1 1
0 = At + (1 − γ)rA + Ay (f + σ 2 )y + Ayy σ 2 y 2
2 2
σ2 2 y A2y (1 − γ)σ 2 yAy (f +
2 σ2
A(1 − γ)(f + 2 − r) 2 − r)(1 − γ)
+ + + . (102)
2γσ 2 2Aγ γ

At this point, it is not obvious how to solve equation (102) due to its non-linearity.
Instead, we will solve it numerically to gain some insights into it.

In our initial experiment, we would like to solve for A(t, y) numerically in the domain
(t, y) ∈ [0, T ] × (0, C) where C is a positive constant. A natural terminal condition to
be applied is
1
A(T, y) = (103)
1−γ
according to the utility function in equation (95). To better analyze the behavior as the
price of risky asset approaches to 0, we apply the change of variable formula x = log(y)
 46
and rewrite equation (102) accordingly as

1 1
0 = At + (1 − γ)rA + Ax (f + σ 2 ) + (Axx − Ax )σ 2
2 2
2 σ2
A(1 − γ)(f + σ2 − r)2 A2x (1 − γ)σ 2 Ax (f + 2 − r)(1 − γ)
+ + + , (104)
2γσ 2 2Aγ γ

where A is understood as A(t, ex ). We then numerically solve equation (104) using an

explicit finite difference scheme. Since our focus here is not the numerical algorithm,
we postpone the detail description of it to later sections. We set T = 0.5 where A is
constant according to equation (103), and iteratively compute the value of A for t < T
and. Plots of A at some earlier times are given in Figure 6 (labeled with different colors).
We observe that for a fixed x, A increases quickly as we stepped back in time. This
can also be visualized by the fact the a curve corresponding to a smaller time is higher
than a curve corresponding to a larger time. At a fixed point in time, represented by
each single curve, A would also increases as x deviates from 0. It is worth mentioning
that when we try to carry out the experiment for larger value of T , numerical blow-up
occurs.

Figure 6: Numerical solution of A at some specific times

Our initial observations motivates us to reduce the scale of function A so that we

can deal with moderate numerical values. To that purpose, we define a new function B
 47
as
B(t, x) = log A(t, ex ), (105)

where x = log y as defined previously. This transformation leads to the a partial differ-
ential equation equivalent to equation (104) in terms of t, x, and B(t, x)

1 1
0 = Bt + (1 − γ)r + Bx (f + σ 2 ) + (Bxx + Bx2 − Bx )σ 2
2 2
2 σ2
(1 − γ)(f + σ2 − r)2 B 2 (1 − γ)σ 2 Bx (f + − r)(1 − γ)
+ + x + 2
, (106)
2γσ 2 2γ γ
µ0 σ 2 +σ02 x
where f = σ 2 +tσ02
after transforming from y to x, together with the terminal condi-
tion B(T, x) = − log(1 − γ). Equation (106) has two advantages over equation (104).
The first is the range of B(t, x) is significantly smaller compare to A(t, y) in equation
(104). The second is the only non-linear term in equation (106) is Bx2 while A(t, y) also
appears in a denominator in equation (104).

A numerical study of equation (106) is conducted using the same algorithm which
was used in the study of equation (104) (the details are postponed to later sections).
Figure 7 displays plots of numerical solution at several point in time. As we can see,
the range of B(t, x) is significantly reduced and we can plot for a much wider range of
x.

Figure 7: Numerical solution of B(t, x) at some specific times

 48
In addition to function B(t, x), we also examined first order numerical derivative
w.r.t. x, which are computed using upwind difference scheme. As shown in Figure 8
(and verified by second order numerical derivatives), the numerically value of Bx (t, x)
changes linearly w.r.t to x. This finding strongly suggests that if we fix the value of t,
then B(t, x) becomes a quadratic function of x. Thus, we are motivated to use a natural
ansatz to solve equation (106) which is illustrated in the next theorem.

Figure 8: Numerical derivative of B(t, x) w.r.t x at some specific times

Theorem 2. Let a, b, and c be functions that map R to R. Assume they solve the
following system of ordinary differential equations on [0, T ] with boundary conditions
a(T ) = 0, b(T ) = 0, and c(T ) = − log(1 − γ). Then B(t, x) = a(t)x2 + b(t)x + c(t) solves
 49
the partial differential equation (106) with terminal condition B(T, x) = − log(1 − γ).

ȧ + 2σ02 2σ 2 2 (1−γ)σ04
2 +tσ 2 ) a + γ a + 2γσ 2 (σ 2 +tσ 2 )2 = 0,



 γ(σ 0 0
σ02 2σ 2 µ0 σ 2 2 −2r 2

+ σ 2γ + r − σ2 )+



 ḃ + [ γ(σ 2 +tσ02 )
+ γ a]b + 2a( γ(σ 2 +tσ02 )

 2 2
 (1−γ)σ ( µ0 σ − r + σ2 ) = 0,
2

0

γσ 2 (σ 2 +tσ02 ) (σ 2 +tσ02 )
(107)
µ0 σ 2 0.5σ 2 −r σ2



 ċ + (1 − γ)r + ( γ(σ 2 +tσ02 )
+ γ + r − 2 )b+

2
 2 2
1−γ µ0 σ
σ 2 a + σ2γ b2 + 2γσ σ 2
2 ( (σ 2 +tσ 2 ) − r + 2 ) = 0,





 0

a(T ) = 0, b(T ) = 0, c(T ) = − log(1 − γ).


Proof. Substitute B(t, x) = a(t)x2 +b(t)x+c(t) into equation (106) to obtain a quadratic
polynomial of x where the coefficients are functions of t. Set the coefficients to be zeros
and match the degree-zero term to the terminal condition of the partial differential
equation: B(T, x) = − log(1 − γ). We then obtain precisely the system of ordinary
differential equations.

The system of equations (107) looks more complicated than it is. Notice only the
first equation is a non-linear equation due to the term a2 . Once the first equation is
solved, the second and third equations can be solved sequentially as first-order linear
ordinary differential equations. It turns out the first equation in (107) is a Riccati
equation which can be solved via standard methods. Under a mild assumption, system
(107) admits a closed-form solution.

Theorem 3. The system of differential equations (107) admits bounded, continuously

differentiable solutions a, b, and c on [0, T ] if σ02 (1 − γ)T < σ 2 γ.

Proof. Once the first equation in (107) is solved, the second and third equations follow
naturally. For the first equation, we can solve it via a standard technique[18] that
transform it into a linear ordinary differential equation. The solution is
−σ02 1 1
a(t) = + , (108)
2σ σ + tσ02
2 2 2σ 2
γ t + CT
where CT is a constant determined by the boundary condition at t = T . In this case,
1−γ σ04
a(T ) = 0 yields CT = 4γσ 4 (σ 2 +T σ02 )2
. Thus, solution to the first equation is

−σ02 1 1
a(t) = + . (109)
2σ 2 σ 2 + tσ02 2σ 2 2σ 4
+ 2σ 2 T (1 − γ1 )
γ t + σ02
 50
While σ 2 + tσ02 is positive for t ≥ 0, the denominator of the second term on the right-
hand-side of equation (109) can become 0 and leads to a blow-up. Denote h(t) =
2σ 2 2σ 4
γ t + σ02 + 2σ 2 T (1 − γ1 ). Since h(t) is monotone increasing in t, it assumes 0 on [0, T ]
if and only if h(0) ≤ 0, which happens exactly when σ02 (1 − γ)T ≥ σ 2 γ. Thus, as
long as σ02 (1 − γ)T < σ 2 γ, there exist a bounded continuously differentiable solution
to a. Otherwise, a blows up to infinity in finite time. In particular, the inequality is
always satisfied when γ > 1, and we only need to worry about this criteria if γ < 1.
Solutions to the second and third follow naturally due to standard linear ODE theory.
The boundedness is a simple consequence of the fact that a, b, and c are continuous
functions on the closed interval [0, T ].

According to theorem 1 and theorem 2, we can derive a solution to the HJB equation,
together with an expression for the optimal control αt∗ . We summarize the result in the
following theorem.

Theorem 4. There is a closed-form solution V (t, y, p) to the HJB equation(97) when

σ02 (1 − γ)T < σ 2 γ,
2 +b(t) log y+c(t)
V (t, y, p) = p1−γ ea(t)(log y) , (110)

where a(t), b(t), and c(t) are described in Theorem 2. Moreover, the optimal choice of
αt in the HJB equation is given by
 2a(t) σ02   r 1 b(t) µ0 
αt∗ = + log y − − − − . (111)
γ γσ 2 σ02 t + γσ 4 γσ 2 2γ γ γσ 2 + γσ02 t

In particular, αt∗ is a linear function of log y for any fixed t.

Proof. Notice B(t, x) and A(t, y) are the same function expressed in different scales.
According to Theorem 2, B(t, x) = a(t)x2 + b(t)x + c(t) solves equation (106). Using
the coordinate transforms x = log y, B(t, x) = log A(t, ex ), and equation (101), we
obtain the expression for V (t, y, p). Then the expression for αt∗ follows from equation
(99).

Now that we have figured out the solution to the HJB equation, what’s left is to
prove the solution is indeed the value function of the original stochastic control problem.
In addition, we also want to verify αt∗ is the optimal control which produced the value
 51
function. For the proof, we will use the HJB equation and adapt a classic Change-of-
Numeraire method[35].

Theorem 5. Among all Gaussian control processes with bounded mean and variance
on the interval [0, T ] that are also jointly Gaussian with with the log-price process Xt ,
and T satisfies the criteria in Theorem 3. Then V (t, y, p) and αt∗ in Theorem 4 are the
value function and the optimal control of the stochastic control problem, respectively.

Proof. Consider at time t ∈ [0, T ], the price of risky asset is y > 0 and the portfolio
value is p > 0. For s ∈ [t, T ], denote the process Yst,y the price process whose value is
y at time t. Also, denote the process Pst,y,p,α the portfolio process where its value is p
at time t, the price of risky asset is y at time t, and the control process α is applied
on the interval [t, T ]. Apply Itô’s lemma to the process V (s, Yst,y , Pst,y,p,α ) and collect
terms with ds and dWs , we get

dV (s, Yst,y , Pst,y,p,α ) = Ls ds + Hs dWs , (112)

where

∂V ∂V 1
Ls = (s, Yst,y , Pst,y,p,α ) + (s, Yst,y , Pst,y,p,α )(f + σ 2 )Yst,y
∂t ∂y 2
∂V h 1 i 1 ∂2V
+ (s, Yst,y , Pst,y,p,α ) αs (f + σ 2 − r) + r Pst,y,p,α + (s, Yst,y , Pst,y,p,α )σ 2 (Yst,y )2
∂p 2 2 ∂y 2
1 ∂2V t,y t,y,p,α 2 2 t,y,p,α 2 ∂2V
+ (s, Y s , Ps )σ α t (P s ) + (s, Yst,y , Pst,y,p,α )αt σ 2 Yst,y Pst,y,p,α ,
2 ∂p2 ∂p∂y
and
∂V ∂V
Hs = (s, Yst,y , Pst,y,p,α )σYst,y + (s, Yst,y , Pst,y,p,α )αs σPst,y,p,α .
∂y ∂p
Notice the terms in Ls are those in the original HJB equation (97) only with stochastic
arguments and we always have Ls ≤ 0 and the equality is achieved if αs = αs∗ .

Next we perform a Change-of-Numeraire using Hs and then apply Girsanov theorem[20].

Using the expression of V (t, y, p) given in Theorem 4, we can deduce
 
Hs = V (s, Yst,y , Pst,y,p,α ) 2a(s) log Yst,y + b(s) + (1 − γ)αs σ. (113)
 52
Since the logarithm of Yst,y and αs are jointly Gaussian with bounded mean and variance
on [0, T ], so is their sum. Set
 
Js = 2a(s) log Yst,y + b(s) + (1 − γ)αs σ. (114)

By Fernique’s theorem[10], there is  > 0 such that

h R (u+)∧T 1 2
i
J dτ
E exp u 2 τ < +∞ (115)

for any u ∈ [0, T ]. By dividing [0, T ] into finite many sub-intervals with length 2 , the
Novikov condition is satisfied according to Corollary 3.5.14 in Karatzas and Shreve[20].
Consequently, the Doleans exponential of Js is a martingale for s ∈ [t, T ]
Rs Rs
Jτ dWτ − 12 Jτ2 dτ
ηst = exp t t . (116)

V (s,Yst,y ,Pst,y,p,α )
By applying Itô’s lemma to ηst (use ηst as the Numeraire)we get

 V (s, Y t,y , P t,y,p,α )  L

s s s
d t
= t ds, (117)
ηs ηs

which is non-positive and is (always) zero when αs∗ is applied. By comparing s = t and
s = T we see
V (t, y, p)ηTt ≥ V (T, YTt,y , PTt,y,p,α ). (118)

Now take expectation on both side, using the fact that V (t, y, p) is non-random, ηTt is
p1−γ
a martingale, and V (T, y, p) = 1−γ we see

h (P t,y,p,α )1−γ i
T
V (t, y, p) ≥ E , (119)
1−γ
where the right hand side is the expected utility function at T .

Notice when α∗ is applied on [t, T ], Ls will always be zero and so is the right-hand-
side of equation (117). In this case, equality is achieved in equation (118). Finally note
αs∗ is a Gaussian process with bounded mean and variance on [0, T ] thanks to the fact
Xst,y = log(Yst,y ) and equation (88). As the deterministic functions in equation (114)
are bounded and continuous, the proof is concluded.
 53
From Theorem 4 and Theorem 5 we see the value function V exist on [0, T ] for a
suitable value T . If T is too large, V will reach infinity in finite time due to the blow-up
of a(t). For a fix time t, V is unbounded and grows quickly as the log-price deviates
from zero. An intuitive explanation consists of two parts. The first being the drift
of price process is determined by current price through equation (87) in the Bayesian
framework. That is, if the price increases then we get a larger drift, and the price is
likely to get even larger due to the increased drift. The second is an investor is allowed
to take arbitrary long/short positions in both assets. We will see that the value function
will behave differently if we restrict the investment action one can take.

6.5 Stochastic Control Problem - Single Action

In this section, we consider a special stochastic control problem with the same market
dynamics described in equations (90) through (94), and the value function defined in
equations (95) and (96). The only difference is that we restrict the the investment ac-
tion to be a single point, i.e. A = {α} for a fixed α ∈ R. Obviously, the HJB equation
becomes simpler as we do not need to maximize a quadratic function within the equa-
tion. Yet this problem provides useful insight regarding the boundary behavior when
the price is very large or close to zero. In later section we will see how it helps us specify
the correct boundary conditions in numerical simulations.

Since there is only one choice of αt ≡ α, the HJB equation (97) is reduced to

1 1
0 = Vt + Vy (f + σ 2 )y + Vp [α(f + σ 2 − r) + r]p
2 2
1 1
+ Vyy σ 2 y 2 + Vpp α2 σ 2 p2 + Vyp ασ 2 yp. (120)
2 2
We assume the value function has the same form as defined in equation (101) (for the
same reason we used to write out equation (101))

V (t, y, p) = A(t, y)p1−γ . (121)

 54
After applying transformations x = log y and B(t, x) = log A(t, ex ), we obtain an equiv-
alent formulation of equation (120) in terms of x and B(t, x)

1 1
0 = Bt + Bx (f + σ 2 ) + σ 2 [Bxx + (Bx )2 − Bx ] + (1 − γ)r
2 2
1 1
+ (1 − γ)(−γ)σ 2 α2 + Bx (1 − γ)σ 2 α + (1 − γ)(f + σ 2 − r)α. (122)
2 2
As in the unconstrained control problem, we conduct a numerical experiment to get
some idea about how does the solution behave. The experiment is based on equation
(122) and explicit finite difference scheme. We fix α = 0.5 in this experiment. The
underlying region is (t, x) ∈ [0, 0.5] × [−100, 100]. The terminal condition applied at
t = 0.5 is B(0.5, x) = − log(1 − γ). Plots of numerically calculated B(t, x) can be found
in Figure 9. It suggests that for a fixed t the function B(t, x) is a linear function of x.

Figure 9: Numerical solution of B(t, x) at some specific times (α = 0.5)

It turns out the numerical study again provides us with useful insight about the value
function in the singleton control problem. We summarize the finding in the following
theorem.

Theorem 6. Let b, and c be functions that map R to R. Assume they solve the following
system of ordinary differential equations on [0, T ] with boundary conditions b(T ) = 0,
and c(T ) = − log(1 − γ). Then B(t, x) = b(t)x + c(t) solves the partial differential
 55
equation (122) with terminal condition B(T, x) = − log(1 − γ).

σ02 σ02
ḃ + σ2 +σ 2 t b + (1 − γ)α σ 2 +σ 2 t = 0,



 0 0

ċ + b[ µ2 0 σ22 + (1 − γ)σ 2 α] + 1 σ 2 b2 + (1 − γ)r+


σ +σ t 2
0 (123)
1 2 2 µ0 σ 2 1 2
2 (γ − 1)γσ α + (1 − γ)( σ 2 +σ02 t + 2 σ − r)α = 0,






b(T ) = 0, c(T ) = − log(1 − γ).


Proof. Substitute B(t, x) = b(t)x + c(t) into equation (122) to obtain a linear function
of x where the coefficients are functions of t. Set the coefficients to be zeros and
match the degree-zero term to the terminal condition of the partial differential equation:
B(T, x) = − log(1−γ). We then obtain the system of ordinary differential equations.

Notice the first equation in equation (123) is a first order linear ordinary differential
equation of b, which can be easily solved. Once we obtain b, c can be solved by inte-
gration according the second equation. Thus, the existence of solution to the system is
established.

Theorem 7. The system of differential equations (123) admits bounded, continuously

differentiable solutions b, and c on [0, T ].

Proof. The existence and uniqueness of solutiona is guaranteed by the Picard–Lindelöf

Theorem[22]. Being continuous functions on an closed interval, the solutions must be
bounded as well.

As a result of Theorem 6 and Theorem 7, we confirm that the function B(t, x) is

indeed a linear function in x for a fixed time t. In addition, we are able to obtain a
closed-form solution to the HJB equation and the solution is the value function in the
singleton control problem.

Theorem 8. For any T ≥ 0, there is a closed-form solution V (t, y, p) to the HJB

equation(36)
V (t, y, p) = p1−γ eb(t) log y+c(t) , (124)

where b(t) and c(t) are as described in theorem 5. In addition, V (t, y, p) is the value
function in the singleton control problem with A = {α}.
 56
Proof. We can follow the same verification argument of Theorem 5, or a direct applica-
tion of Feynman-Kac formula[33] will lead to the conclusion.

We conclude this section by discussing how the value function behave and compare
it to the value function in the unconstrained control problem. To make the discussion
easier, we factor out p1−γ and work in the (t, x, B(t, x)) system.

The first thing we can do is to go back and check the accuracy of the numerical
solution illustrated in Figure 9. It turns out the numerical solution provides good ap-
proximations to the close-form solution. In particular, for a fixed time t, B(t, x) is
linearly increasing with the log-price of the risky asset. This can be partly explained
by the rule that an investor must always puts a fixed portion of wealth in the risky
asset. Since the drift of risky asset is positively related to its price, the expected value
drops/increases when the price drops/increases. In comparison, for the unconstrained
control problem, investors can expect large value function no matter high or low the
observed prices are, since they are allowed to take arbitrary positions in both assets.

The second concerns the limit behavior of B(t, x) when x approaches +∞ or −∞.
From Theorem 2, we see function B(t, x) grows quadratically in x for any fixed time t
in the unconstrained control problem (Figure 7). In addition, the optimal control αt∗
is a linear function of x for a fixed t, describing how an investor adjusts the position
in risky asset according to the underlying price. In the singleton control problem, an
invest always put a fixed portion of wealth in the risky asset and the function B(t, x) is
a linear function of x. Of interest is a the case in-between: a control problem where an
investor has some but limited flexibility. For example, we can set the action set A to be
a closed interval [l, h]. Intuitively, an investor can adjust the optimal decision αt∗ within
the interval based on Xt , the log-price of risky asset. For a large positive Xt , αt∗ will be
close to or even equal to h since the risky asset is expected to give a good return. For a
large negative Xt , αt∗ is likely to be close or equal to l for a similar reason. Once h (or l)
has been reached, even one observes larger (or smaller) values of Xt , no further actions
are allowed to be taken due to the constraint A = [l, h]. Thus, for extreme values of x,
the investor is essentially taking a constant position (l or h) in the risky asset. Based
on our study of the singleton control problem. We can expect the value functions to be
 57
asymptotically linear in x for a fixed time t.

Remark. In the case where A is a set of finite cardinality or closed interval, we do not
expect to be able to derive a closed-form solution. Instead, we will rely on numerical
study to obtain a numerically smooth solution. From this perspective, previous discussion
helps us determine the proper boundary conditions.

6.6 Stochastic Control Problem - Finite Action Set

Based on the discussion on the single action problem, we can make a generalization by
allowing the action set A to contain finite number of values. The main purpose is to
build a connection to the constrained action problems (section 6.7) where the action set
is a closed interval. Intuitively, when we approximate an closed interval with a (fine)
finite set, the resulting constrained stochastic control problem should have similar solu-
tion to its finite action counterpart. We will see that is indeed the case based on result
of numerical study, but in this section we focus on the finite action problem.

When the action set is finite, an investor can choose how much wealth to allocate in
the risky asset but the choices are limited and discrete (compare to the unconstrained
problem). For example, if we set A = {0, 1}, then an investor can choose to invest only
in bond or only in stock. Of course A can contain values greater than 1 or smaller than
0, representing potential short positions in two assets. As we shall see, the solutions are
not fundamentally different from the single action case.

Assume the action set contains k distinct real numbers, i.e. A = {ai |1 ≤ i ≤ k, ai ∈
R}. Since the actions set is finite, the HJB equation (97) is equivalent to
n 1 1
0 = Vt + max Vy (f + σ 2 )y + Vp [αt (f + σ 2 − r) + r]p
αt ∈A 2 2
1 1 o
+ Vyy σ 2 y 2 + Vpp αt2 σ 2 p2 + Vyp αt σ 2 yp . (125)
2 2
On the right hand side of equation (125), there is a discrete optimization problem over
a finite set. We do not expect to obtain a closed-form solution to the the optimal choice
 58
αt . However, in numerical study it can be solved by simply comparing among all pos-
sible values. To reduce the dimension and scale of the problem, apply transformations
V (t, y, p) = A(t, y)p1−γ , x = log y, and B(t, x) = log A(t, ex ), we reach the following
equation

1 1
0 = Bt + Bx (f + σ 2 ) + σ 2 [Bxx + (Bx )2 − Bx ] + (1 − γ)r+
2 2
n 1 1 o
(1 − γ)max − γσ 2 αt2 + (Bx σ 2 + f + σ 2 − r)αt . (126)
αt ∈A 2 2
The optimal choice αt∗ is determined by comparison among all possible values:
 
∗ 1 2 2 2 1 2
αt = argmax − γσ α + (Bx σ + f + σ − r)α . (127)
α∈A 2 2

Equation (127) is useful when it come to solve equation (126) numerically since it re-
duces an optimization problem into finding the maximal element in an array, which can
be integrated into classic finite difference schemes.

Although the nature of the discrete optimization problem makes it difficult to find
a closed-form solution to equation (126), we can still pursue numerical solutions thanks
to the finiteness of the action set A. The idea is to apply simple explicit finite difference
methods while solving optimization problems via equation (127) at each grid point. Be-
fore proceeding to numerical study, we discuss the terminal and boundary conditions.
It is obvious that at time T , B(T, x) = − log(1 − γ). The boundary conditions as
x → ±∞ deserve some discussion. As seen in the single action problem, B(t, x) is a
linear function of x if t is fixed. In addition, it is intuitive to allocation more (less) in
the risky asset if x is large (small) due to the market dynamics. Thus, in finite action
problems, we expected the optimal control to be constant for extreme values of x, and
the function B(t, x) is asymptotically linear in x in the limit (which agrees with the
result obtained without boundary conditions).

We outline the numerical algorithm before presenting result of numerical studies.

Start by evenly discretizing the time domain [0, T ] and price domain [−X, X] into M +1
and 2N + 1 grid points, respectively. Then we assign terminal condition at time T by
setting B(T, x) = − log(1 − γ). After that, an iteration is carried out to compute
 59
B(t, x), using explicit finite difference scheme based on equation (126). In addition, the
optimal control αt∗ at each grid point is computed via equation (127).For grid points
corresponding to x = X or x = −X, a linear boundary condition is applied to determine
B(t, x). The complete process is summarized in the algorithm below.
Algorithm 2: Finite Difference Scheme for Finite Action Problems
Input: T , X, M , N . A = {0, 1}, γ = 0.5, µ0 = 0.03, r = 0.03, σ = 0.3, and
σ0 = 0.03.
1 Discretize [0, T ] and [−X, X] into M and 2N sub-intervals.
(i) Denote ti = iT /M for 0 ≤ i ≤ M .
(ii) Denote xj = (j − N )X/N for 0 ≤ j ≤ 2N .
2 Assign terminal condition:
B(tM , xj ) = − log(1 − γ) for 0 ≤ j ≤ 2N .
3 For each j between M − 1 and 0:
(i) For each i between 1 and 2N − 1:
(i-a) Compute Bx (tj , xi ) and Bxx (tj , xi ) using B(tj+1 , xi−1 ), B(tj+1 , xi ), and
B(tj+1 , xi+1 ) via finite difference.
(i-b) Compute αt∗ via equation (127).
(i-c) Compute B(tj , xi ) via equation (126).
(ii) Apply linear boundary conditions:
B(tj , x0 ) = 2B(tj , x1 ) − B(tj , x2 ), and
B(tj , x2N ) = 2B(tj , x2N −1 ) − B(tj , x2N −2 ).
We present some case studies of the finite action control problem. For the purpose
of comparison, we restrict numerical studies in (t, x) ∈ [0, 2] × [−1, 1]. We start with the
simplest case where A = {0, 1}, meaning an investor can invest only in the risky asset
or only in the risk-free asset at a given time. Plots of numerical solutions to B(t, x)
at different times are shown in Figure 10, together with the optimal actions. In Figure
10(a) we see plots of B(t, x) tend to be flat for small (negative) values of x, and the
corresponding optimal controls are αt∗ = 0. This phenomenon agrees with earlier dis-
cussion which suggests that when the price of risky asset is small it tends to be even
smaller in the near future. Consequently, it is safer to invest all the wealth in bond
rather than stock. In addition, B(t, x) behaves like its linear counterpart in the single
action problem where an investor invest in the bond all the time. On the other hand, as
 60
x increases, B(t, x) becomes upward-sloping and the optimal controls switch to 1, which
agrees with our discussion for exactly the same reason. We also observe that B(t, x)
tend to be linear when x approaches to 1.
 61

(a) Numerical solution of B(t, x)

(b) Optimal control

Figure 10: Numerical solution of B(t, x) and optimal control

 62
We can examine the asymptotic behavior of B(t, x) in more detail by comparing plots
of numerical solution with the plots of analytic solutions where A = {0} and A = {1}.
In Figure 11, color plots represent the numerical solution we obtained from Algorithm
2, where A = {0, 1}. The smaller t is, the higher the curve corresponding to t lies
in the figure. The dashed black lines represent analytical solutions obtained following
Theorem 6. More specifically, flat dashed black lines represent analytic solutions where
A = {0}, and only left halves are shown to make the comparison more clear. Dashed
black lines with positive slopes represent analytic solutions where A = {1}, and only
right halves are displayed. We can clearly see in Figure 11 that the numerical solutions
almost coincide with analytic solutions when x becomes very large or small. In addition,
as x increase from −1 to 1, numerical solutions change gradually from horizontal lines
on the left to upward-sloping lines on the right.

Figure 11: Comparison of the numerical solution with analytical solutions to single
action problems (A = {0} and A = {1})
 63
In addition to asymptotic properties, we are also interested in the smoothness of
the numerical solutions. Based on numerical solution to B(t, x), we can compute its
first-order and second-order derivatives w.r.t x. In our studies, numerical derivatives are
computed using central finite difference schemes. As for boundary points, the deriva-
tives are implied by linear boundary condition. Figure 12(a) shows plots of first-order
numerical derivatives at different times. It suggests that Bx (t, x) gradually increases
from 0 which represents all-bond investment, to different positive levels that represent
all-stock investment. The further a moment is from terminal time, the higher the level
on the right-hand-side is. This agrees with earlier result in the section of single action
problems. Figure 12(b) shows plots of second-order numerical derivatives at correspond-
ing times. We see that the slope of B(t, x) is always increasing but the rate it changes
is not a constant. Starting from x = −1, the first-order derivative increases at a grow-
ing rate until a point near zero. After that the rate it increases keeps declining. In
summary, our analysis shows how the numerical solutions to B(t, x) are connected to
the solutions to the single action problems, and suggests they are numerically twice
continuously differentiable w.r.t x, which is a property implicitly assumed by the HJB
equation (97).
 64

(a) 1st order derivative Bx (t, x)

(b) 2nd order derivative Bxx (t, x)

Figure 12: Numerical derivatives of B(t, x)

 65
Next we add some flexibility to the previous finite action problem by setting A =
{−1, 0, 1}. In this case, an investor can choose to take a short position in the risky asset
with αt = −1. Consequently, the expected return can be substantial even for smaller
(negative) values of x by simply shorting the risky asset and invest proceedings in the
risk-free asset. Notice with this modification, equation (126), equation (127), and the
numerical procedure outlined in Algorithm 2 remain valid. In fact, the only difference
in numerical study is that the computer needs to make extra comparisons due to the
increased size of A.

We present the result of numerical studies with the same setup except for the admis-
sible action set A. First we examine the numerical solutions to B(t, x) and corresponding
optimal actions. Numerical solutions to B(t, x) at different times are plotted in Figure
13(a). The right half of Figure 13(a) is similar to that of Figure 8(a), but the left part
is quite different. Since an investor can choose to short the stock when x is small, the
left halves of plots are no longer flat. We can get a better idea about this from Figure
13(b) which shows optimal actions one would take at different times w.r.t. different x.
For small or large x, e.g. x < −0.2 or x > 0.1, an investor would take a −100% short
position or a 100% long position in the stock. Note for a small interval in between -0.2
and 0.1, one would put everything in the bond, i.e. choosing αt∗ = 0. The transition is
by no means smooth since the admissible action set itself is not continuous.
 66

(a) Numerical solution of B(t, x)

(b) Optimal control

Figure 13: Numerical solution of B(t, x) and optimal control

 67
Next we compare the numerical solution to analytic solutions of single action prob-
lems where the control is set to be A = {−1} and A = {1} , respectively. As the optimal
control takes either the smallest or largest possible value in A when x is small enough or
large enough, we expect the numerical solution to B(t, x) to be asymptotically similar to
solutions of those single action control problems. We can visualize that effect in Figure
14. Similar to Figure 11, numerically calculated solutions are plotted in solid color lines
and analytic solutions are plotted in dashed black lines. For large or small x, curves
with color are not distinguishable with dashed line segments which represent analytical
solutions where A = {−1} or A = {1}, respectively.

Figure 14: Comparison of the numerical solution with analytical solutions to single
action problems (A = {−1} and A = {1})

Numerical derivatives w.r.t x are displayed in Figure 15. Figure 15(a) shows how do
the slopes of B(t, x) gradually change from lower levels corresponding to single action
problem with A = {−1} to higher levels corresponding to single action problem with
 68
A = {1}. Figure 15(b) shows the corresponding second order derivatives. Notice there
are relatively plat regions on top of curves( comparing to those in Figure 12(b). Overall,
the transitional pattern w.r.t x is similar to what we have seen in Figure 7. The result
again suggests the numerical solution is well-behaved w.r.t. x.
 69

(a) 1st order derivative Bx (t, x)

(b) 2nd order derivative Bxx (t, x)

Figure 15: Numerical derivatives of B(t, x)

 70
We conclude this section with a case where the admissible control set is further
refined. Define A = {0.01i − 1|0 ≤ i ≤ 200, i ∈ N} and it’s a superset of {−1, 0, 1}
we used in the previous case study. We want to check what happens when a much
“larger” control set is put into use. In addition, the refined set can be considered as
an discrete approximation to the closed interval [−1, 1]. We hope to connect it with
the constrained stochastic control problems in the next section. Without reiterating the
numerical process, we present results in Figure 16, 17, and 18, respectively. From Figure
16(a) and Figure 17, we observe that the numerical solution to B(t, x) is very similar
to what displayed in Figure 13(a) and Figure 14. However, optimal control plotted in
Figure 16(b) is distinct from its counterpart in Figure 13(b). As an investor now has
many more choices with the refined control set, the optimal control changes gradually
from -1 to 1 in a linear fashion. The derivative plots in Figure 18 are also similar to
those in Figure 15. The only visible difference is perhaps in the second order numerical
derivatives. Notice in Figure 18(b), the plots tend to be smoother “bell-shaped” curves,
and there is no plat regions on top as in Figure 15.
 71

(a) Numerical solution of B(t, x)

(b) Optimal control

Figure 16: Numerical solution of B(t, x) and optimal control

 72

Figure 17: Comparison of the numerical solution with analytical solutions to single
action problems (A = {−1} and A = {1})
 73

(a) 1st order derivative Bx (t, x)

(b) 2nd order derivative Bxx (t, x)

Figure 18: Numerical derivatives of B(t, x)

 74
6.7 Stochastic Control Problem - Constrained Action
In this section, we consider a stochastic control problem with the same market dynamics
and utility function as defined in previous sections. The only difference is the control
set becomes a closed interval, i.e. A = [m, n] ⊂ R where m and n are constant real
numbers. In particular, when m < 0 or n > 1, an investor can take a limited short
position in the risky or risk-free asset.

Substitute A = [m, n] into equation (97) yields

n 1 1
0 = Vt + sup Vy (f + σ 2 )y + Vp [αt (f + σ 2 − r) + r]p
αt ∈[m,n] 2 2
1 1 o
+ Vyy σ 2 y 2 + Vpp αt2 σ 2 p2 + Vyp αt σ 2 yp . (128)
2 2
Notice the second term on the right-hand-side involves a constrained quadratic opti-
mization problem, and we do not expect to find a closed-form solution for V (t, y, p).
To reduce the number of variables and make it easier to carry out numerical ex-
periments, we again apply transformations V (t, y, p) = A(t, y)p1−γ , x = log y, and
B(t, x) = log A(t, ex ). Then an equivalent formulation in terms of x, t, and B(t, x) is

1 1
0 = Bt + Bx (f + σ 2 ) + σ 2 [Bxx + (Bx )2 − Bx ] + (1 − γ)r+
2 2
n 1 1 o
(1 − γ) sup − γσ 2 αt2 + (Bx σ 2 + f + σ 2 − r)αt . (129)
αt ∈[m,n] 2 2

The quadratic function of αt in the curly bracket attains its maximum when αt =
Bx σ 2 +f +0.5σ 2 −r
γσ 2
. But due to the constraint αt ∈ [m, n], the optimal αt∗ should be deter-
mined as
 Bx σ 2 + f + 0.5σ 2 − r

αt∗ = min n, max m, . (130)
γσ 2
Thus, we can set up numerical experiments similar to the finite-action cases based on
equation (130). More specifically, at any point (t0 , x0 ) ∈ [0, T ]×R, αt∗ can be determined
by the first order numerical derivative w.r.t x, plus comparisons with m and n. In terms
of boundary conditions, since the admissible control set is bounded, asymptotically the
optimal control should be either m or n. Thus, B(t, x) should still be asymptotically
linear in x for any fixed time. The rest of the numerical algorithm is the same as the
 75
finite-action case and we outline the algorithm below.
Algorithm 3: Finite Difference Scheme for Constrained Action Problems
Input: T , X, M , N . A = {0, 1}, γ = 0.5, µ0 = 0.03, r = 0.03, σ = 0.3, and
σ0 = 0.03.
1 Discretize [0, T ] and [−X, X] into M and 2N sub-intervals.
(i) Denote ti = iT /M for 0 ≤ i ≤ M .
(ii) Denote xj = (j − N )X/N for 0 ≤ j ≤ 2N .
2 Assign terminal condition:
B(tM , xj ) = − log(1 − γ) for 0 ≤ j ≤ 2N .
3 For each j between M − 1 and 0:
(i) For each i between 1 and 2N − 1:
(i-a) Compute Bx (tj , xi ) and Bxx (tj , xi ) using B(tj+1 , xi−1 ), B(tj+1 , xi ), and
B(tj+1 , xi+1 ).
(i-b) Compute αt∗ via equation (130).
(i-c) Compute B(tj , xi ) via equation (129).
(ii) Apply linear boundary conditions:
B(tj , x0 ) = 2B(tj , x1 ) − B(tj , x2 ), and
B(tj , x2N ) = 2B(tj , x2N −1 ) − B(tj , x2N −2 ).
We present some result of numerical studies based on Algorithm 3. For the purpose
of comparison, we restrict numerical studies to (t, x) ∈ [0, 2] × [−1, 1], which is the same
as the finite-action problems.

First we consider the case where A = [0, 1]. Plots of numerical solutions to B(t, x)
and optimal actions are shown in Figure 19. In Figure 19(a), we see the plots gradually
transit from flat to positive-slope. The corresponding optimal controls are mostly 0 on
the left part and 1 on the right part. However, in the middle part of Figure 19(b) we
observe optimal controls transit from 0 to 1. Although the interval where the changes
take place is small, it makes a fundamental difference from what we observe in Figure
10(b) where optimal controls jump from 0 to 1.
 76

(a) Numerical solution of B(t, x)

(b) Optimal control

Figure 19: Numerical solution of B(t, x) and optimal control

 77
In addition, we plot numerical solutions against analytic solutions to single-action
problems where the actions are held constant, i.e. A = {0} or A = {1}. In Figure 20,
color plots represent our numerical solutions and dashed black plots represent analytic
solutions where A = {0} and A = {1}, respectively. Clearly they match very well for
large or small x.

Figure 20: Comparison of the numerical solution with analytical solutions to single
action problems (A = {0} and A = {1})

We can further examine numerical derivatives of numerical solution w.r.t. x. As

shown in Figure 21(a), the first order derivatives gradually transit from 0 to a constant
values, which are slopes of solutions to the single-action problems. From Figure 21(b)
we see the second order derivatives are bell-shaped curves which are similar to Figure
18(b) in the finite-action section.
 78

(a) 1st order derivative Bx (t, x)

(b) 2nd order derivative Bxx (t, x)

Figure 21: Numerical derivatives of B(t, x)

 79
Next we look at the numerical experiment where the admissible control set A =
[−1, 1]. In this case, an investor is allowed to take not only any long-only positions, but
also any short position in the risky asset up to -1. The results are plotted in Figure 22,
23, and 24. Since short positions are allow, B(t, x) is no longer flat when x approaches
to -1, as shown in Figure 22(a). Together with Figure 23, we see the numerical solutions
asymptotically match analytical solutions to constant control problem where the action
is always -1 or 1. In terms of optimal controls, Figure 22(b) shows they change gradually
from -1 to 1 in an interval between -0.2 and 0.2, and assume either -1 or 1 otherwise. In
Figure 24, we see first order numerical derivatives transit smoothly from a set of negative
values to a set of positive value, corresponding to two constant control problems.
 80

(a) Numerical solution of B(t, x)

(b) Optimal control

Figure 22: Numerical solution of B(t, x) and optimal control

 81

Figure 23: Comparison of the numerical solution with analytical solutions to single
action problems (A = {−1} and A = {1})
 82

(a) 1st order derivative Bx (t, x)

(b) 2nd order derivative Bxx (t, x)

Figure 24: Numerical derivatives of B(t, x)

 83
We end this section by discussing how the the optimal control changes over time. In
Figure 22(b) (or similarly Figure 19(b)), the curves representing the change from −1 to
1 (or from 0 to 1) are getting steeper as we step back in time, and the intervals where
the change actually happen is getting smaller. An investor would choose either -1 or 1
outside those shrinking intervals. This effect can be better visualized by looking at the
color map in Figure 25 which corresponding to the case where A = [−1, 1]. In Figure
25, the horizontal axis is the log-price of risky asset, the vertical axis is time, and the
color of each pixel represents the optimal control corresponding to a specific price and a
specific time. The pure blue region on the left means an investor will choose α∗ = −1,
the pure yellow region on the right means an investor will choose α∗ = 1. Between the
pure blue and pure yellow region, the optimal control changes from −1 to 1. As we can
see, the transitional region shrinks when we move from t = 2 to t = 0. The observation
agrees with equation (130) where the fractional term is a linear function of Bx and x
with positive coefficients for a fixed t, and the coefficients become larger if t is smaller
(see equation (91) for detail). Thus, if we fix x and move back in time, it is easier for
the optimal control to hit upper or lower bound since due to the expanding value of Bx
(e.g. Figure 18(a)).
 84

Figure 25: Evolution of optimal control in time (A = [−1, 1])

Will the transitional interval shrink to a single point? We don’t have a complete
answer to that. However, transitional interval can become very small in some cases.
For example, when we conduct a study on (t, x) ∈ [0, 2] × [−5, 5] with A = [−30, 30]
and other parameters the same as previous studies. As we can see in Figure 26, the
transitional region is approximately [−4, 4] near t = 2, but it almost shrinks to a single
point when t is close to 0. As a matter of fact, we discretize the price range [−5, 5] into
501 evenly separated grid points, and there are only four consecutive points (out of 500)
where the optimal control is neither −30 nor 30 when t < 0.1.
 85

Figure 26: Evolution of optimal control in time (A = [−30, 30])

6.8 Discussion
In this chapter, we have studied 4 different cases of the portfolio allocation problem
based on one common stochastic differential equation with Bayesian update. The main
difference among these cases is the action set A. For the unconstrained-action and
single-action problems, we were able to establish theoretical results and find closed-
form solutions. We notice that a problem similar to the unconstrained-action problem
has been studied and one can refer to Bismuth[5] for a more general framework in the
multi-dimensional case. The single-action problem provides us with important infor-
mation about the asymptotic behavior of the value function when the price is extreme,
which gives us the correct boundary conditions in the numerical studies of the finite-
action and constrained-action problems. We then took a deep dive into the finite-action
and constrained-action problems with several interesting cases of action set, and obtain
 86
some interesting numerical results.

In section 6.6, we have conducted experiment for a finite-action problem where the
admissible control set is A = {0.01i − 1|0 ≤ i ≤ 200, i ∈ N}. The results are presented
in Figure 16, 17, and 18. When compare those results to what we obtained from the
constrained-action problem using A = [−1, 1], we do not see much difference. In other
words, numerical solutions in two cases are very close and share the same asymptotic
behavior. In fact, the maximal difference in B(t, x) is negligible based on our computa-
tion. In addition, optimal controls and numerical derivatives also behave very similarly
in those two cases. We did notice that the actions in the finite-action problem must
be integer multiples of 0.01 and there is no such constraint in the constrained-action
problem. Despite the difference, the finite-action problem is qualified to be a very close
approximation to the constrained-action problem according to our result.

To conclude this chapter, we compare the unconstrained-action problem with the

classic Merton’s portfolio problem. According to Merton, if the drift of risky asset µ is
known, an investor should always hold a constant portion of risky asset (equation (62)).
In our model, we do not know the exact value of µ. Under Bayesian belief system, we
are able to reach equation equation (87) which suggests that the drift of the risky asset
changes dynamically with its price. If we substitute the drift term of equation (87) into
equation (62), we get the “optimal control”

µ0 σ 2 + σ02 x r
ω∗ = 2 2 2
− 2 , (131)
(σ + tσ0 )σ γ σ γ

where x stands for the realization of Xt at time t. Using above equation and equation
(111), we can directly compute

2a(t) b(t) 1
α∗ − ω ∗ = x+ + , (132)
γ γ 2γ
where the relation x = log y is used. Thus, we see that applying Merton’s classic result
with dynamic drift only gives sub-optimal investment choice. In particular, α∗ and ω ∗
differ by a linear function of x for a fixed t.
Chapter 7

Conclusion

We have studied a variety of portfolio optimization models in this thesis. Started with
the classic Markowitz model, we went through several modifications related to robust
optimization and Bayesian methods. Then we took a dynamic perspective and examined
portfolio optimization models in multiple periods and continuous-time. By combining a
continuous-time model with Bayesian update, we derived a stochastic differential equa-
tion whose drift changes dynamically according observed price. Based on the equation,
we further studied different cases of optimal portfolio allocation problems, and obtain
theoretical as well as numerical results.

There are two directions worth further exploration. The first is to study the exis-
tence of classic solutions to HJB equations in finite-action and constrained-action cases.
Although we have obtained numerically smooth solutions in both cases, the study is
not completely rigorous due to the lack of theoretical support. The second direction
is to apply effective numerical methods for HJB equations. We have relied on special
coordinate transformations that make HJB equations easier to be studied, with linear or
quadratic boundary conditions. However, such transformation may be hard to find for
more general stochastic control problems and utility functions. As existing numerical
methods often requires strong regularity conditions and limited boundary conditions ,
a general numerical framework about HJB equation will be much appreciated.

87
References

[1] Anderson, T. W. An introduction to multivariate statistical analysis. John Wiley

& Sons, New York, 1984.

[2] Bain, A., and Crisan, D. Fundamentals of stochastic filtering, vol. 60. Springer
Science & Business Media, 2008.

[3] Bellman, R. The theory of dynamic programming. Bulletin of the American

Mathematical Society 60, 6 (1954), 503–516.

[4] Best, M. J., and Grauer, R. R. On the sensitivity of mean-variance-efficient

portfolios to changes in asset means: some analytical and computational results.
The review of financial studies 4, 2 (1991), 315–342.

[5] Bismuth, A., Guéant, O., and Pu, J. Portfolio choice, portfolio liquidation, and
portfolio transition under drift uncertainty. Mathematics and Financial Economics
13, 4 (2019), 661–719.

[6] Black, F., and Litterman, R. Global portfolio optimization. Financial Ana-
lysts Journal 48 (1992), 28–43.

[7] Black, F., and Litterman, R. B. Asset allocation: Combining investor views
with market equilibrium. The Journal of Fixed Income 1, 2 (1991), 7–18.

[8] Bolstad, W. M., and Curran, J. M. Introduction to Bayesian statistics. John

Wiley & Sons, 2016.

[9] Broadie, M. Computing efficient frontiers using estimated parameters. Annals

of operations research 45, 1 (1993), 21–58.

88
 89
[10] Da Prato, G., and Zabczyk, J. Stochastic equations in infinite dimensions.
Cambridge university press, 2014.

[11] Fama, E. F., and French, K. R. Common risk factors in the returns on stocks
and bonds. Journal of Financial Economics 33 (1993), 3–56.

[12] Fama, E. F., and French, K. R. A five-factor asset pricing model. Journal of
financial economics 116, 1 (2015), 1–22.

[13] Folland, G. B. Real analysis: modern techniques and their applications, vol. 40.
John Wiley & Sons, 1999.

[14] Friedman, A. Stochastic differential equations and applications. In Stochastic

differential equations. Springer, 2010, pp. 75–148.

[15] Gilbarg, D., and Trudinger, N. S. Elliptic partial differential equations of

second order, vol. 224. springer, 2015.

[16] Goldfarb, D., and Iyengar, G. Robust portfolio selection problems. Mathe-
matics of operations research 28, 1 (2003), 1–38.

[17] Grant, M., and Boyd, S. Cvx: Matlab software for disciplined convex program-
ming, version 2.1.

[18] Hille, E. Ordinary differential equations in the complex domain. Courier Corpo-
ration, 1997.

[19] Iyengar, G. N. Robust dynamic programming. Mathematics of Operations Re-

search 30, 2 (2005), 257–280.

[20] Karatzas, I., and Shreve, S. Brownian Motion and Stochastic Calculus,
vol. 113. Springer Science & Business Media, 2012.

[21] Keeping, E. S. Introduction to statistical inference. Courier Corporation, 1995.

[22] Knapp, A. W. Basic real analysis. Springer Science & Business Media, 2005.

[23] Krylov, N. V. Nonlinear elliptic and parabolic equations of the second order,
vol. 7. Springer, 1987.
 90
[24] Lee, W. Theory and methodology of tactical asset allocation, vol. 65. John Wiley
& Sons, 2000.

[25] Lintner, J. The valuation of risk assets and the selection of risky investments
in stock portfolios and capital budgets. The Review of Economics and Statistics
(1965), 13–37.

[26] Lu, Z. Robust portfolio selection based on a joint ellipsoidal uncertainty set.
Optimization Methods & Software 26, 1 (2011), 89–104.

[27] Markowitz, H. Portfolio selection. The Journal of Finance, 07 (1952), 77–91.

[28] Merton, R. C. Lifetime portfolio selection under uncertainty: The continuous-

time case. The review of Economics and Statistics (1969), 247–257.

[29] Merton, R. C. An analytic derivation of the efficient portfolio frontier. Journal

of financial and quantitative analysis (1972), 1851–1872.

[30] Merton, R. C., et al. Optimum consumption and portfolio rules in a continuous-
time model. Journal of Economic Theory 3, 4 (1971), 373–413.

[31] Meucci, A. Risk and Asset Allocation. Springer Science & Business Media, 2007.

[32] Meucci, A. The black-litterman approach: Original model and extensions.

Shorter version in, THE ENCYCLOPEDIA OF QUANTITATIVE FINANCE,
Wiley (2010).

[33] Pham, H. Continuous-time stochastic control and optimization with financial ap-
plications, vol. 61. Springer Science & Business Media, 2009.

[34] Sharpe, W. F. Capital asset prices: A theory of market equilibrium under con-
ditions of risk. The journal of finance 19, 3 (1964), 425–442.

[35] Shreve, S. E. Stochastic calculus for finance II: Continuous-time models, vol. 11.
Springer Science & Business Media, 2004.

[36] Treynor, J. L. Market value, time, and risk. Time, and Risk (August 8, 1961)
(1961).
 91
[37] Treynor, J. L. Toward a theory of market value of risky assets. 1961.
Appendix A

A.1 Derivation of HJB Equation

We follow the notations used in Chapter 3 of Pham[33]. Consider a n−dimensional
stochastic process Xs ∈ Rn defined by

dXs = b(s, Xs , αs )ds + σ(s, Xs , αs )dWs , (133)

where s stands for time, αs is the control process, and Ws is a d−dimensional Brownian
motion. At some future time T ∈ R, a reward will be given according to a function g(·)
with XT as the input to g(·). Assume at time t < T the observed value of Xs is x ∈ Rn ,
we can define the value function of (t, x) as

v(t, x) = sup E[g(XTt,x )]. (134)

α∈A

In equation (134), XTt,x stands for the stochastic process whose value is x at time t,
with the control α applied between t and T , and A is the set of all admissible control
processes. The objective of a stochastic control problem is the find the optimal control
process α∗ so that the expectation on the right-hand-side of equation (134) is maximized.

Now consider applying a constant control αs = a between t and t + h where a

is an proper choice w.r.t. the pair (t, x). By definition of v(t, x) and the dynamic
programming principle (section 3.3 of Pham[33]), we have the following relation

t,x
v(t, x) ≥ E[v(t + h, Xt+h )]. (135)

92
 93
Assume v is smooth enough, we may apply Ito’s formula between t and t + h:
Z t+h
t,x ∂v
v(t + h, Xt+h ) = v(t, x) + ( + La v)(s, Xst,x )ds + Ms , (136)
t ∂t
where
1
La v = b(t, x, a) · Dx v + tr(σ(t, x, a)σ(t, x, a)T Dxx v) (137)
2
and Ms is a (local) martingale. In equation (137), Dx and Dxx stand for the first order
and second order derivative operator, respectively. In general, we do not need to worry
about Ms since eventually we will let h →
− 0 and the expectation of Ms can be handled
by localization. Thus, substituting equation (136) into equation (135) we have
Z t+h 
∂v a t,x
0≥E ( + L v)(s, Xs )ds . (138)
t ∂t

Divide both sides of equation (138) by h and let h →

− 0, we obtain

∂v
0≥ (t, x) + La v(t, x). (139)
∂t
Now since a can be any admissible choice, we have
∂v
0≥ (t, x) + sup [La v(t, x)] , (140)
∂t a∈A

where A represents the set of admissible choices of a w.r.t. (t, x). However, if a is chosen
according to the optimal control that achieves equation (134), using the same argument
we obtain the equality
∂v
0= (t, x) + sup [La v(t, x)] (141)
∂t a∈A

assuming the supreme is finite. Finally, since (t, x) is arbitrary, we obtain the HJB
equation
∂v
0= (t, x) + sup [La v(t, x)] . (142)
∂t α∈A