BLS WORKING PAPERS

U.S. DEPARTMENT OF LABOR

Bureau of Labor Statistics

 

OFFICE OF PRICES AND LIVING

CONDITIONS

Higher Moments in Perturbation Solution of the Linear-Quadratic

Exponential Gaussian Optimal Control Problem

Baoline Chen, Rutgers University

Peter A. Zadrozny, U.S. Bureau of Labor Statistics

Working Paper 348

November 2001

 

The views expressed are those of the authors and do not necessarily reflect the policies of the U.S. Bureau of Labor

Statistics or the views of other staff members. We thank David Belsley and Larry Karp for comments.
HIGHER MOMENTS IN PERTURBATION SOLUTION OF THE LINEAR-QUADRATIC EXPONENTIAL

GAUSSIAN OPTIMAL CONTROL PROBLEM*

Baoline Chen

Department of Economics

Rutgers University

Camden, NJ 08102

tel: (856)-225-6025

e-mail: baoline@crab.rutgers.edu

Peter A. Zadrozny

Bureau of Labor Statistics

2 Massachusetts Ave., NE, Room 3105

Washington, DC 20212

tel: (202)-691-6591

e-mail: zadrozny_p@bls.gov

March 21, 2001

ABSTRACT

The paper obtains two principal results. First, using a new definition of

higher-order (>2) matrix derivatives, the paper derives a recursion for

computing any Gaussian multivariate moment. Second, the paper uses this result

in a perturbation method to derive equations for computing the 4th-order Taylorseries approximation of the objective function of the linear-quadratic

exponential Gaussian (LQEG) optimal control problem. Previously, Karp (1985)

formulated the 4th multivariate Gaussian moment in terms of MacRae’s definition

of a matrix derivative. His approach extends with difficulty to any higher (>4)

multivariate Gaussian moment. The present recursion straightforwardly computes

any multivariate Gaussian moment. Karp used his formulation of the Gaussian 4th

moment to compute a 2nd-order approximation of the finite-horizon LQEG objective

function. Using the simpler formulation, the present paper applies the

perturbation method to derive equations for computing a 4th-order approximation

of the infinite-horizon LQEG objective function. By illustrating a convenient

definition of matrix derivatives in the numerical solution of the LQEG problem

with the perturbation method, the paper contributes to the computational

economist’s toolbox for solving stochastic nonlinear dynamic optimization

problems.

*The paper represents the authors’ views and does not represent any official

positions of the U.S. Bureau of Labor Statistics. We thank David Belsley and

Larry Karp for comments. Forthcoming in Computational Economics.
1

 1. INTRODUCTION.

Consider the discrete-time state equation and feedback control rule

(1.1) xt = F0xt-1 + G0ut + εt,

(1.2) ut = P0xt-1,

where x is an n×1 state vector, u is an m×1 control variable, and ε is an n×1

disturbance ~ NIID(0,Σ), F0 and G0 are n×n and n×m parameter matrices, and P0 is

an m×n feedback control matrix. The discounted linear-quadratic (LQ) objective

function is

(1.3) v(xt-1,N) = (1/2)∑ =

N

k 0 δk

( T xt+k−1 Q0xt+k-1 + 2 T ut+k S0xt+k-1 + T ut+k R0ut+k),

where 0 < δ ≤ 1 is the discount factor and N is a finite or infinite horizon,

Q0, S0, and R0 are given n×n, m×n, and m×m preference parameter matrices. We

assume that the quadratic form in (1.3) is non-negative definite overall and

positive definite with respect to u. Superscript T denotes vector or matrix

transposition.

The risk-avoiding linear-quadratic-exponential-Gaussian (LQEG) objective

function is

(1.4) J(xt-1,N) = Et-1exp[-v(xt-1,N)],

where Et-1 denotes expectation conditioned on variables realized in period t-1.

The discrete-time LQEG problem is: given xt-1, N, and the parameters F0, G0, Σ,

R0, S0, and Q0, minimize (1.4) with respect to P0, subject to (1.1) to (1.3).

Jacobson (1973) showed that, for finite N, the optimal P0 is obtained by

iterating on a discrete-time recursive Riccati equation. As N → ∞, the equation

converges to a nonrecursive or algebraic Riccati equation, which can be solved

quickly and accurately using the Schur-decomposition method (Laub, 1979). See

also Hansen and Sargent (1995). Karp (1985) addressed the problem of determining

the contribution of higher moments (> 2) of ε to the value of the LQEG objective

function, J(xt-1,N). Expanding J(xt-1,N) in a Taylor-series and using some matrix

differentiation rules of MacRae (1974), Karp derived an algorithm for computing
2

a two-term approximation of J(xt-1,N) based on the 2nd and 4th moments of ε (odd

moments of ε are zero). The complexity of MacRae’s differentiation rules

apparently dissuaded Karp from attempting to derive equations for computing

higher-order (> 2) terms based on higher moments (> 4) of ε.

The present paper extends Karp’s results in two ways. First, using a

simpler representation of matrix derivatives based on the total-differential

rather than on partial-derivative forms (Magnus and Neudecker, 1988), the paper

derives a simple recursion for computing any moments of a Gaussian random

vector. Second, using this result, the paper applies the perturbation method

(Judd, 1998, chs. 13-14) to derive equations for computing the 4th-order Taylorseries approximation of J(xt-1,∞), based on the 2nd and 4th moments of ε. The

second result illustrates using higher Gaussian moments in a perturbation

solution of a nonlinear dynamic stochastic model. There is a growing interest in

economics in solving dynamic stochastic models with the perturbation method

(Anderson and Hansen, 1996; Chen and Zadrozny, 2000b; Collard and Juillard,

2000; Sims, 2000).

The paper proceeds as follows. Section 2 states and proves a theorem and

corollary that give the recursion for computing any multivariate Gaussian

moment. Section 3 derives linear, perturbation-solution equations for computing

the 4th-order approximation of J(xt-1,∞) for the undiscounted problem. Section

4 illustrates the results of Section 3 numerically. Section 5 gives concluding

remarks. There are two technical appendices. Appendix A explains the definitions

and rules of matrix differentiation that are used in Sections 2 and 3. Appendix

B explains how to compute the 1st- to 4th-order derivatives of f(x) =

exp[-(1/2)xT

Qx], in the gradient forms defined in appendix A, that are inputs in

the perturbation solution equations of Section 3.

2. RECURSION FOR COMPUTING ANY MULTIVARIATE GAUSSIAN MOMENT.

This section follows the definitions and rules of matrix differentiation

explained in appendix A, which the reader should read before proceeding.

Let x be a random n-vector distributed N(μ,Σ) and let m(z) = E[exp(xT

z)],

for z ∈ Rn

, be its moment generating function given by

(2.1) m(z) = exp[μT

z + (1/2)zTΣz].
3

As explained in appendix A, ∇k

m(z) is the kth-order gradient matrix of kth-order

partial derivatives of m(z). Let μk be the kth uncentered moment of x defined as

μk = E(Πk⊗x), where Πk⊗x denotes k-1 successive Kronecker products of x (e.g.,

Π2⊗x = x⊗x). Then, μk = vec[∇k

m(0)]. The following theorem states a recursion

for computing ∇k

m(z) for any finite k.

THEOREM 1: Suppose x ∈ Rn

 ~ N(μ,Σ), with moment generating function m(z) given

by (2.1). Then, for k = 3, ..., K,

(2.2) vec[∇k

m(z)] = vec{[∇k-1m(z)T ⊗ (μ + Σz)] + (k-1)[∇k-2m(z)T ⊗ Σ)]},

where ∇k

m(z) is nk-1×n, starting with ∇m(z) = m(z)(μT + zTΣ) and ∇2m(z) = (μ +

Σz)∇m(z) + m(z)Σ.

Theorem 1 immediately implies

COROLLARY 1: For k = 3, ..., K,

(2.3) μk = vec{[mat(μk-1)T ⊗ μ] + (k-1)[mat(μk-2)T ⊗ Σ]},

starting with μ1 = μ and μ2 = (μ ⊗ μ) + vec(Σ), where mat(μk) denotes the

unvectorization of μk to an nk-1×n dimensional matrix.

There is a subtlety in the role of the vectorization operator in (2.3).

One might think we could unvectorize (2.3) and write it as mat(μk) = [mat(μk-1)T

⊗ μ] + (k-1)[mat(μk-2)T ⊗ Σ], but this cannot be done because, whereas mat(μk) is

nk-1×n, [mat(μk-1)T ⊗ μ] + (k-1)[mat(μk-2)T ⊗ Σ] is n2

×nk-2.

Corollary 1 is the principal result of this section. For k = 4, (2.3)

corresponds to Karp’s (1985) equation (10), based on MacRae’s (1974) definition

of a matrix derivative. Whereas (2.3) is valid for any k, it would be very

tedious to extend Karp’s equation (10) correspondingly. We now prove theorem 1.

PROOF OF THEOREM 1:
4

We repeatedly vectorize and differentiate expressions by applying rules

(A.3), (A.19), and (A.20), without referencing them explicitly. We also

repeatedly use the fact that m(z) is a scalar.

1. Derivation of ∇m(z): Differentiating (2.1), we obtain

(2.4) dm(z) = m(z)(μT

 + zTΣ)dz.

Then, because dm(z) = ∇m(z)dz holds for all dz, we can drop dz and obtain ∇m(z)

= m(z)(μT

 + zTΣ).

2. Derivation of ∇2

m(z): Differentiating (2.4) and vectorizing terms, we

obtain

(2.5) d2

m(z) = dm(z)(μT + zTΣ)dz + m(z)dzTΣdz

 = dzT

[(μ + Σz)dm(z) + m(z)Σdz]

 = dzT

[(μ + Σz)∇m(z) + m(z)Σ]dz.

Then, because d2

m(z) = dzT∇m(z)dz holds for all dz, we can drop dz and obtain

∇2

m(z) = (μ + Σz)∇m(z) + m(z)Σ.

3. Derivation of (2.2) for k = 3: Differentiating (2.5) and vectorizing

terms, we obtain

(2.6) d3

m(z) = d2m(z)(μT + zTΣ)dz + 2dm(z)dzTΣdz

 = dzT

[(μ + Σz)d2m(z) + 2Σdzdm(z)]

 = dzT

[(μ + Σz)dzT∇2m(z) + 2Σdz∇m(z)]dz

 = (Π2⊗dzT

)vec[(μ + Σz)dzT∇2m(z) + 2Σdz∇m(z)]

 = (Π2⊗dzT

){[∇2m(z)T

⊗ (μ + Σz)] + 2[∇m(z)T⊗ Σ]}dz.

Then, because d3

m(z) = (Π2⊗dzT)∇3m(z)dz holds for all dz, we can drop dz and

obtain ∇3

m(z) = [∇2m(z)T

⊗ (μ + Σz)] + 2[∇m(z)T⊗ Σ].

4. Derivation of (2.2) for k = 4: Differentiating (2.6) and vectorizing

terms, we obtain
5

(2.7) d4

m(z) = d3m(z)(μT + zTΣ)dz + 3d2m(z)dzTΣdz

 = dzT

[(μ + Σz)d3m(z) + 3Σdzd2m(z)]

 = dzT

[(μ + Σz)(Π2⊗dzT)∇3m(z) + 3ΣdzdzT∇2m(z)]dz

 = (Π2⊗dzT

)vec[(μ + Σz)(Π2⊗dzT)∇3m(z) + 3ΣdzdzT∇2m(z)]

 = (Π2⊗dzT

){[∇3m(z)T ⊗ (μ + Σz)] + 3[∇2m(z)T

⊗ Σ]}(Π2⊗dz)

 = (Π4⊗dzT

)vec{[∇3m(z)T ⊗(μ + Σz)] + 3[∇2m(z)T

⊗ Σ]}.

Then, because d4

m(z) = (Π3⊗dzT)∇3m(z)dz = (Π4⊗dzT)vec[∇3m(z)] holds for all dz,

we can drop dz and obtain vec[∇4

m(z)] = vec{[∇3m(z)T

⊗ (μ + Σz)] + 3[∇2

m(z)T

⊗

Σ]}.

5. Derivation of (2.2) for k ≥ 5: Continuing in this fashion,

(2.8) dk

m(z) = dk-1m(z)(μT + zTΣ)dz + (k-1)dk-2m(z)dzTΣdz

 = dzT

[(μ + Σz)dk-1m(z) + (k-1)Σdzdk-2m(z)]

 = dzT

[(μ + Σz)(Πk-2⊗dzT)∇k-1m(z) + (k-1)Σdz(Πk-3⊗dzT)∇k-2m(z)]dz

 = (Π2⊗dzT

)vec[(μ + Σz)(Πk-2⊗dzT)∇k-1m(z)+ (k-1)Σdz(Πk-3⊗dzT)∇k-2m(z)]

 = (Π2⊗dzT

){[∇k-1m(z)T ⊗ (μ + Σz)] + (k-1)[∇k-2m(z)T

⊗ Σ]}(Πk-2⊗dz)

 = (Πk⊗dzT

)vec{[∇k-1m(z)T ⊗ (μ + Σz)] + (k-1)[∇k-2m(z)T

⊗ Σ]}.

Then, because dk

m(z) = (Πk⊗dzT)vec[∇km(z)] holds for all dz, we can drop dz and

obtain vec[∇k

m(z)] = vec{[∇k-1m(z)T

⊗ (μ + Σz)] + (k-1)[∇k-2m(z)T⊗ Σ]}.

Thus, we have proved theorem 1.

The particular value here of the corollary is that it expresses Gaussian

moments in the same gradient form of matrix derivatives that is used in the

perturbation solution equations in Section 3. From an analytic standpoint, any

arrangement of elements of matrix derivatives is equally acceptable. However, in
6

order to derive useful matrix-algebraic solution equations, the disturbance

moments must be expressed compatibly with matrix derivatives in the perturbation

equations.

3. PERTURBATION SOLUTION OF THE INFINITE-HORIZON LQEG PROBLEM.

We now focus exclusively on the N = ∞ problem and express (1.1) to (1.4)

more compactly.

Using the feedback-control rule (1.2) to eliminate the control vector u

and restating the disturbance as θεt, where θ ≥ 0 is a scalar parameter, we

write the state law of motion, (1.1), in the closed-loop form

(3.1) xt = Φxt-1 + θεt,

where Φ = F0 + G0P0 is the closed-loop matrix. The perturbation method produces

an approximate Taylor-series solution of the stochastic LQEG problem

(henceforth, "stochastic solution"), centered on the exact solution of the

nonstochastic LQEG problem (henceforth, "nonstochastic solution"). When θ = 0,

the LQEG problem is nonstochastic. In the perturbation method, as θ goes from

zero to a positive value it "extrapolates" the nonstochastic solution to an

approximate stochastic solution. When θ > 0 but ≠ 1, the original

parameterization is maintained by rescaling Σ as Σ/θ.

For simplicity, we drop the time subscript. Specifically, we write xt-1,

εt, and Et-1 as x, ε, and E. Then, for N = ∞, we express (1.4) in the manner of

the Bellman equation of dynamic progamming, as

(3.2) g(x) = f(x)Eg(Φx + θε)δ

,

where f(x) = exp[-(1/2)xT

Qx] and Q = P0

T

R0P0 + P0

T

S0 + S0

T

P0 + Q0, so that g(x) =

J(x,∞). For given values of x and the parameters, hence, for given f(x), Φ, θ,

and Σ/θ, the g(x) function which solves (3.2) is the stochastic solution. Hansen

and Sargent (1995) consider a related logarithmic form of (3.2).

When θ = 0, (3.2) reduces to

(3.3) g0(x) = f(x)g0(Φx).
7

Whereas the analytical expression of the stochastic solution is unknown, the

nonstochastic solution, g0(x), which solves (3.3) for given values of x and the

parameters, is given by

(3.4) g0(x) = exp[-(1/2)xT

Vx],

where V = ∑∞

k =0 δk

(ΦT)kQΦk.

We assumed that the quadratic form in (1.3) is non-negative definite

overall and positive definite with respect to u, which implies that Q is

positive definite (Q > 0). We now also assume that Φ is a stable matrix, i.e.,

has all eigenvalues less than one in absolute value. The latter assumption is

concisely stated as ρ(Φ) < 1, where ρ(Φ) denotes the spectral radius of Φ,

i.e., the largest absolute eigenvalue of Φ.

Given Q > 0 and 0 < δ ≤ 1, ρ(Φ) < 1 is sufficient for V to exist (be

finite) and Q > 0 implies V > 0. By writing V = Q + δΦT

[Q + δΦTQΦ + δ

2

(ΦT)2QΦ2 +

...]Φ and noting that the expression in brackets is the same as V, we see that

V = ΦT

VΦ + Q. Then, vectorizing this equation using rule (A.3), we obtain

(3.5) vec(V) = [ 2 nI – δ(Π2⊗ΦT

)]-1vec(Q),

where 2 nI denotes the n2

×n2 identity matrix. Because 0 < δ ≤ 1, the assumption

ρ(Φ) < 1 implies δρ(Π2⊗Φ) < 1, hence, that 2 nI – δ(Π2⊗ΦT) is nonsingular.

Our objective is to compute the 4th-order approximate solution of the

infinite-horizon problem without discounting (δ = 1),

(3.6) gˆ(x) = go + ∇go(x–xo) + (1/2)(x–xo)T

∇2go(x–xo)

 + (1/6)[Π2⊗(x–xo)T

]∇3go(x-xo) + (1/24)[Π3⊗(x–xo)T]∇4go(x-xo),

where x0 is any value of the state vector and g0 = g(x0), ..., ∇4

g0 = ∇4g(x0). We

restrict the discussion to the undiscounted problem for the sake of brevity. The

discounted problem is treated in exactly the same way but involves considerably

more algebra.
8

Because the perturbation method treats θ as a state variable, we define

the current and next period’s augmented state vectors as z = (xT

,θ)T and z′ =

xT

ΦT + θεT, θ)T, and write (3.2) as

(3.7) g*

(z) = f(x)Eg*(z′)δ.

The 4th-order approximation of g*

(z) is

(3.8) gˆ (z) * = *g0 + * ∇g0 (z–zo) + (1/2)(z–zo)T *

0

2 ∇ g (z-zo)

 + (1/6)[Π2⊗(z–zo)T

] *

0

3 ∇ g (z-zo) + (1/24)[Π3⊗(z–zo)T

] *

0

4 ∇ g (z-zo),

where z0 = (x0

T

,0)T and *g0 = g*(z0), ..., *

0

4 ∇ g = ∇4

g*(z0).

The perturbation method proceeds as follows: (1) we set δ = 1 and

differentiate (3.7) four times in succession with respect to z; (2) we convert

the differential form equations into gradient form; (3) we match up coefficients

of like powers of elements of dz, which indicates solution equations for *g0 ,

..., *

0

4 ∇ g ; (4) we set θ = 1 in (3.8), and consolidate terms into the form of

(3.6).

In the following four subsections, we skip most steps in the derivations

because including them all would make the paper unbearably long. Full notes on

the derivations are available from the authors. However, most steps are

elementary, direct or inverse, applications of the vectorization rule (A.3). The

direct applications of (A.3) are based on the observations that a scalar is an

example of a column vector and that a column vector equals the vectorization of

itself. The steps rely mostly on rules (A.3), (A.19), and (A.20), and are used

without explicit referencing. The solution equations for ∇ *g0 , ..., *

0

4 ∇ g derived

in this section require ∇f(x), ..., ∇f4

(x) as input values. The equations for

computing these gradients are derived in appendix B. We use the fact that odd

moments of the normally distributed disturbance vector are zero. Finally, for

brevity, we write the undiscounted (3.7) as g*

 = fEg*′, and its differentials

correspondingly.

3.1. FIRST-ORDER PERTURBATION-SOLUTION EQUATIONS.
9

Setting δ = 1 and differentiating (3.7) with respect to z, we obtain

(3.9) dg*

 = dfEg*′ + fEdg*′.

Corresponding to the partition dz = (dxT

, dθ)T, the columns of the

gradient of g*

 partition as ∇g* = [ * ∇g1 , * ∇g2 ]. Then, writing the differentials in

terms of their gradients and evaluating expectations, we obtain

(3.10) * ∇g1 dx + * ∇g2 dθ = [g*

′∇f + f * ∇g1 Φ]dx + f * ∇g2 dθ.

Because (3.10) holds for all dx and dθ, we can drop dx and dθ and obtain

(3.11) * ∇g1 = * ∇g1 fΦ + g*

′∇f,

(3.12) * ∇g2 = f * ∇g2 .

Note that the * ∇g1 ’s on the left and right sides of (3.11) are identical

and evaluated at z0 = (x0

T

, 0)T, even though they, respectively, stem from g* and

g*

′. The reason is that applying the chain rule of differentiation and

evaluating the expectation in (3.9) implements the prime in dg*

′, so that it can

be omitted from (3.10) to (3.12) and all derivatives of g*

 are evaluated at z0.

The same thing happens in the higher-order perturbation equations. By contrast,

the prime stays on g*

′, which is evaluated at z0 as g0(Φx0) = exp[-(1/2)x0

T

ΦTVΦx0].

The maintained assumptions, Q > 0 and ρ(Φ) < 1, imply that f < 1 and

fρ(Φ)k

 < 1, hence, that knI - f(Πk⊗ΦT) is invertible, for k = 0, 1, 2, ... .

Generalizing this section’s results, we would see that the invertibility of kn

I

- f(Πk⊗ΦT), for k = 0, 1, ..., K, is sufficient and necessary for the existence

and uniqueness of the coefficients *g0 , ..., *

0

K ∇ g of the Kth-order approximate

solution.

At this point, given that In – fΦ is invertible, (3.11) and (3.12) imply

(3.13) * ∇g1 = g*

′∇f[In - fΦ]-1,
10

and * ∇g2 = 01×1.

3.2. SECOND-ORDER PERTURBATION-SOLUTION EQUATIONS.

Differentiating (3.9) with respect to z, we obtain

(3.14) d2

g* = d2fEg*′ + 2dfEdg*′ + fEd2g*′.

Corresponding to the partition dz = (dxT

, dθ)T, the rows and columns of

the 2nd-order gradient of g*

 partition as

(3.15) ∇2

g* = 















∇ ∇

∇ ∇

*

22

* 2

21

2

*

12

* 2

11

2

g g

g g

.

Theoretically, ∇2

g* is symmetric because cross-partial derivatives of g*

are equal. However, practically, we require only that (z–zo)T∇2g*(z–zo) is

correct regardless of the structure of ∇2

g*. That is, a quadratic form can be

defined variously, in terms of a symmetric matrix, an upper or lower triangular

matrix, or one with no particular pattern. The perturbation method produces n

linear restrictions for determining the 2n elements of *

21

2 ∇ g and *

12

2 ∇ g , so that n

additional linear restrictions must be introduced to determine *

21

2 ∇ g and *

12

2 ∇ g

uniquely. The simplest restrictions are *

21

2 ∇ g = 01×n or *

12

2 ∇ g = 0n×1. Arbitrarily,

we choose *

21

2 ∇ g = 01×n.

Then, writing the differentials in (3.14) in terms of their gradients,

imposing the previous zero solution ( * ∇g2 = 01×1) and the current zero

restriction ( *

21

2 ∇ g = 01×n), vectorizing terms, and evaluating expectations, we

obtain

(3.16) (Π2⊗dxT

)vec( *

11

2 ∇ g ) + dxT

dθ *

12

2 ∇ g + (dθ)2 *

22

2 ∇ g

 = (Π2⊗dxT)[vec(g*′∇2f + 2∇fT * ∇g1 Φ) + f(Π2⊗ΦT)vec( *

11

2 ∇ g )]

 + dxT

dθfΦT *

12

2 ∇ g + (dθ)2

[fvec(Σ)Tvec( *

11

2 ∇ g ) + f *

22

2 ∇ g ],
11

where the first term in the second brackets on the right side of the equation is

obtained from EεT *

11

2 ∇ g ε = Evec(εT ⊗ εT

)vec( *

11

2 ∇ g ) = vec(Σ)T

vec( *

11

2 ∇ g ).

Because (3.16) holds for all dx and dθ, we can drop dx and dθ and obtain

(3.17) vec( *

11

2 ∇ g ) = f(Π2⊗ΦT

)vec( *

11

2 ∇ g ) + vec(g*

′∇2f + 2∇fT * ∇g1 Φ),

(3.18) *

12

2 ∇ g = fΦT *

12

2 ∇ g ,

(3.19) *

22

2 ∇ g = f *

22

2 ∇ g + fvec(Σ)T

vec( *

11

2 ∇ g ).

Then, given that kn

I - f(Πk⊗ΦT) is invertible for k = 0, 1, 2, (3.17) to

(3.19) imply

(3.20) vec( *

11

2 ∇ g ) = [ 2n

I – f(Π2⊗ΦT)]-1vec(g*′∇2f + 2∇fT * ∇g1 Φ),

(3.21) *

22

2 ∇ g = f(1 – f)-1vec(Σ)T

vec( *

11

2 ∇ g ),

and *

12

2 ∇ g = 0n×1.

3.3. THIRD-ORDER PERTURBATION-SOLUTION EQUATIONS.

Differentiating (3.14) with respect to z, we obtain

(3.22) d3

g* = d3fEg*′ + 3d2fEdg*′ + 3dfEd2g*′ + fEd3g*′.

Corresponding to the partitions Π2⊗zT

 = (Π2⊗dxT, dxTdθ, dxTdθ, (dθ)2) and

dz = (dxT

, dθ)T, the rows and columns of the 3rd-order gradient of g* partition

as
12

(3.23) ∇3

g* = 









































∇ ∇

∇ ∇

∇ ∇

∇ ∇

*

42

* 3

41

3

*

32

* 3

31

3

*

22

* 3

21

3

*

12

* 3

11

3

g g

g g

g g

g g

.

Considerations of redundancy that led us to set *

21

2 ∇ g = 01×n now lead us to set

*

21

3 ∇ g = *

31

3 ∇ g = 0n×n, *

41

3 ∇ g = 01×n, and *

32

3 ∇ g = 0n×1.

Then, writing the differentials in (3.22) in terms of their gradients,

imposing previous and current zero restrictions and zero solutions, vectorizing

terms, and evaluating expectations, we obtain

(3.24) (Π3⊗dxT

)vec( *

11

3 ∇ g ) + (Π2⊗dxT

)dθ *

12

3 ∇ g + dxT

(dθ)2 *

22

3 ∇ g + (dθ)3 *

42

3 ∇ g

 = (Π3⊗dxT){f(Π3⊗ΦT)vec( *

11

3 ∇ g ) + vec[g*

′∇3f + 3(ΦT * ∇g1

T ⊗ ∇2

f)

 + 3(ΦT *

11

2 ∇ g Φ ⊗ ∇fT

)]} + (Π2⊗dxT)dθf(Π2⊗ΦT) *

12

3 ∇ g

 + dxT(dθ)2{fΦT *

22

3 ∇ g + 3∇fT

vec(Σ)Tvec( *

11

2 ∇ g ) + 3∇fT *

22

2 ∇ g

 + fΦT *

11

3 ∇ g T

vec(Σ) + f[vec(Σ)T ⊗ ΦT]( 2 n,n K + 3n I )vec( *

11

3 ∇ g )}

 + (dθ)3

[f *

42

3 ∇ g + fvec(Σ)T *

12

3 ∇ g ].

The last term in (3.24), involving the permutation matrix, 2 n,n K , is derived

using (A.5).

Because (3.24) holds for all dx and dθ, we can drop dx and dθ and obtain

(3.25) vec( *

11

3 ∇ g ) = f(Π3⊗ΦT

)vec( *

11

3 ∇ g ) + vec[g*

′∇3f + 3(ΦT * ∇g1

T ⊗ ∇2

f)

 + 3(ΦT *

11

2 ∇ g Φ ⊗ ∇fT

)],
13

(3.26) *

12

3 ∇ g = f(Π2⊗ΦT

) *

12

3 ∇ g ,

(3.27) *

22

3 ∇ g = fΦT *

22

3 ∇ g + 3∇fT

vec(Σ)Tvec( *

11

2 ∇ g ) + 3∇fT *

22

2 ∇ g + fΦT *

11

3 ∇ g T

vec(Σ)

 + f[vec(Σ)T ⊗ ΦT

]( 2 n,n K + 3n I )vec( *

11

3 ∇ g ),

(3.28) *

42

3 ∇ g = f *

42

3 ∇ g + fvec(Σ)T *

12

3 ∇ g .

Then, given that kn

I - f(Πk⊗ΦT) is invertible for k = 0, 1, 2, 3, (3.25)

to (3.28) imply

(3.29) vec( *

11

3 ∇ g ) = [ 3n

I – f(Π3⊗ΦT)]-1vec[g*′∇3f + 3(ΦT * ∇g1

T ⊗ ∇2

f)

 + 3(ΦT *

11

2 ∇ g Φ ⊗ ∇fT

)],

(3.30) *

22

3 ∇ g = [In - fΦT

]-1{3∇fTvec(Σ)Tvec( *

11

2 ∇ g ) + 3∇fT *

22

2 ∇ g + fΦT *

11

3 ∇ g T

vec(Σ)

 + f[vec(Σ)T ⊗ ΦT

]( 2 n,n K + 3n I )vec( *

11

3 ∇ g )},

*

12

3 ∇ g = 1 2n 0 × , and ∇3

g42 = 01×1.

3.4. FOURTH-ORDER SOLUTION EQUATIONS.

Differentiating (3.22) with respect to z, we obtain

(3.31) d4

g* = d4fEg*′ + 4d3fEdg*′ + 6d2fEd2g*′ + 4dfEd3g*′ + fEd4g*′.

Corresponding to the partitions Π3⊗zT

 = (Π3⊗dxT, (Π2⊗dxT)dθ, (Π2⊗dxT)dθ,

dxT

(dθ)2, (Π2⊗dxT)dθ, dxT(dθ)2, dxT(dθ)2, (dθ)3) and dz = (dxT, dθ)T, the rows and

columns of the 4th-order gradient of g*

 partition as

(3.32) ∇4

g*T =
14

 





















∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇

∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇

* T

82

* T 4

72

* T 4

62

* T 4

52

* T 4

42

* T 4

32

* T 4

22

* T 4

12

4

* T

81

* T 4

71

* T 4

61

* T 4

51

* T 4

41

* T 4

31

* T 4

21

* T 4

11

4

g g g g g g g g

g g g g g g g g

.

As before, considerations of redundancy lead us to set *

21

4 ∇ g = *

31

4 ∇ g = *

51

4 ∇ g =

n 2n 0 × , *

41

4 ∇ g = *

61

4 ∇ g = *

71

4 ∇ g = 0n×n, *

32

4 ∇ g = *

52

4 ∇ g = 1 2n 0 × , *

62

4 ∇ g = *

72

4 ∇ g = 0n×1, and

*

81

4 ∇ g = 01×n.

Then, writing the differentials in (3.31) in terms of their gradients,

imposing previous and current zero restrictions and zero solutions, vectorizing

terms, and evaluating expectations, we obtain

(3.33) (Π4⊗dxT

)vec( *

11

4 ∇ g ) + (Π3⊗dxT

)dθ *

12

4 ∇ g + (Π2⊗dxT

)(dθ)2 *

22

4 ∇ g

 + dxT(dθ)3 *

42

4 ∇ g + (dθ)4 *

82

4 ∇ g = (Π4⊗dxT

){f(Π4⊗ΦT)vec( *

11

4 ∇ g )

 + vec[g*

′∇4f + 4(ΦT * ∇g1

T ⊗ ∇3

f) + 6(ΦT *

11

2 ∇ g Φ ⊗ ∇2

fT)

 + 4(∇fT ⊗ vec((Π2⊗ΦT

)∇3g11Φ))]} + (Π3⊗dxT)dθf(Π3⊗ΦT) *

12

4 ∇ g

 + (Π2⊗dxT)(dθ)2{f(Π2⊗ΦT) *

22

4 ∇ g + 6[vec(Σ)T

vec( *

11

2 ∇ g ) + *

22

2 ∇ g ]vec(∇2

f)

 + 4[(vec(Σ)T ⊗ ΦT

)( 2 n,n K + 3n I )vec( *

11

3 ∇ g ) ⊗ ∇fT

]

 + 4[ΦT *

11

3 ∇ g T

vec(Σ) ⊗ ∇fT] + 4[ΦT *

22

3 ∇ g ⊗ ∇fT

]

 + f[(vec(Σ)T ⊗ (Π2ΦT

))( 3 n,n K + 2 ,n 2n K + 4n I )vec( *

11

4 ∇ g )]

 + f[(ΦT ⊗ vec(Σ)T ⊗ ΦT

)((Kn,n ⊗ 2nI )+( 2nI ⊗ Kn,n)+ 4n I )vec( *

11

4 ∇ g )]}

 + dxT

(dθ)3{fΦT∇4g42 + f(vec(Σ)T ⊗ ΦT)(2 2 n,n K + 3n I ) *

12

4 ∇ g }

 + (dθ)4

{f *

82

4 ∇ g + 3fvec(Σ ⊗ Σ)T

vec( *

11

4 ∇ g ) + fvec(Σ)T *

22

4 ∇ g }.
15

In (3.33), terms involving the permutation matrices, K, are derived using (A.5).

Because (3.33) holds for all dx and dθ, we can drop dx and dθ and obtain

(3.34) vec( *

11

4 ∇ g ) = f(Π4⊗ΦT

)vec( *

11

4 ∇ g ) + vec{g*

′∇4f + 4(ΦT * ∇g1

T ⊗ ∇3

f)

 + 6(ΦT *

11

2 ∇ g Φ ⊗ ∇2

fT) + 4[∇fT ⊗ vec((Π2⊗ΦT)∇3g11Φ)]},

(3.35) *

12

4 ∇ g = f(Π3⊗ΦT

) *

12

4 ∇ g ,

(3.36) *

22

4 ∇ g = f(Π2⊗ΦT

) *

22

4 ∇ g + 6[vec(Σ)T

vec( *

11

2 ∇ g ) + *

22

2 ∇ g ]vec(∇2

f)

 + 4[(vec(Σ)T ⊗ ΦT

)( 2 n,n K + 3n I )vec( *

11

3 ∇ g ) ⊗ ∇fT

]

 + 4[ΦT *

11

3 ∇ g T

vec(Σ) ⊗ ∇fT] + 4[ΦT *

22

3 ∇ g ⊗ ∇fT

]

 + f[(vec(Σ)T ⊗ (Π2ΦT

))( 3 n,n K + 2 ,n 2n K + 4n I )vec( *

11

4 ∇ g )],

 + f[(ΦT ⊗ vec(Σ)T ⊗ ΦT

)((Kn,n ⊗ 2n I )+( 2nI ⊗ Kn,n)+ 4n I )vec( *

11

4 ∇ g )]

(3.37) *

42

4 ∇ g = fΦT

∇4g42 + f(vec(Σ)T ⊗ ΦT)(2 2 n,n K + 3n I ) *

12

4 ∇ g ,

(3.38) *

82

4 ∇ g = f *

82

4 ∇ g + 3fvec(Σ ⊗ Σ)T

vec( *

11

4 ∇ g ) + fvec(Σ)T *

22

4 ∇ g .

Then, given that kn

I - f(Πk⊗ΦT) is invertible for k = 0, 1, 2, 3, 4,

(3.34) to (3.38) imply

(3.39) vec( *

11

4 ∇ g ) = [ 4n I - f(Π4⊗ΦT

)]

-1

vec{g*′∇4f + 4(ΦT * ∇g1

T ⊗ ∇3

f)

 + 6(ΦT *

11

2 ∇ g Φ ⊗ ∇2

fT) + 4[∇fT ⊗ vec((Π2⊗ΦT)∇3g11Φ)]},
16

(3.40) *

22

4 ∇ g = [ 2n I - f(Π2⊗ΦT

)]

-1{6[vec(Σ)T

vec( *

11

2 ∇ g ) + *

22

2 ∇ g ]vec(∇2

f)

 + 4[(vec(Σ)T ⊗ ΦT

)( 2 n,n K + 3n I )vec( *

11

3 ∇ g ) ⊗ ∇fT

]

 + 4[ΦT *

11

3 ∇ g T

vec(Σ) ⊗ ∇fT] + 4[ΦT *

22

3 ∇ g ⊗ ∇fT

]

 + f[(vec(Σ)T ⊗ (Π2ΦT

))( 3 n,n K + 2 ,n 2n K + 4n I )vec( *

11

4 ∇ g )],

 + f[(ΦT ⊗ vec(Σ)T ⊗ ΦT

)((Kn,n ⊗ 2n I )+( 2nI ⊗ Kn,n)+ 4n I )vec( *

11

4 ∇ g )]}

(3.41) *

82

4 ∇ g = f(1 – f)-1[3vec(Σ ⊗ Σ)T

vec( *

11

4 ∇ g ) + vec(Σ)T *

22

4 ∇ g ],

*

12

4 ∇ g = 1 3n 0 × , and *

42

4 ∇ g = 0n×1.

3.5 SUMMARY OF SOLUTION EQUATIONS.

After imposing zero restrictions and zero solutions, (3.8) reduces to

(3.42) gˆ*

(x,θ) = * g0 + * ∇g1 (x–xo) + (1/2){(x–xo)T *

11

2 ∇ g (x-xo) + *

22

2 ∇ g θ2

}

 + (1/6){[Π2⊗(x–xo)T

] *

11

3 ∇ g (x-xo) + (x-x0)T *

22

3 ∇ g θ2

}

 + (1/24){[Π3⊗(x–xo)T

] *

11

4 ∇ g (x-xo) + [Π2⊗(x–xo)T

] *

22

4 ∇ g θ2

 + *

82

4 ∇ g θ4

}.

Setting θ = 1 and consolidating terms in (3.42), writing the result in the

form of (3.6), and matching coefficients with those in (3.6), we obtain

(3.43) go = * g0 + (1/2) *

22

2 ∇ g + (1/24) *

82

4 ∇ g ,

(3.44) ∇go = * ∇g1 + (1/6) *

22

3 ∇ g T

,

(3.45) ∇2

go = *

11

2 ∇ g + (1/24)mat( *

22

4 ∇ g ),
17

(3.46) ∇3

go = *

11

3 ∇ g ,

(3.47) ∇4

go = *

11

4 ∇ g ,

where mat(∇4 * g22 ) denotes the unvectorization of ∇4 * g22 to an n×n matrix, * g0 =

exp[-(1/2)x0

T

Vx0] and * ∇g1 , *

11

2 ∇ g , *

22

2 ∇ g , *

11

3 ∇ g , *

22

3 ∇ g , *

11

4 ∇ g , *

22

4 ∇ g , and *

82

4 ∇ g are

given by (3.13), (3.20), (3.21), (3.29), (3.30), (3.39), (3.40), and (3.41).

Extending the solution to the 5th-order adds the equation ∇5

go = *

11

5 ∇ g and

the terms (1/120) *

16,1

5 ∇ g and (1/120) *

22

5 ∇ g to the right sides of (3.44) and (3.46).

In other words, extending a solution’s order not only adds new solution

equations for higher-order terms but also adds terms to previous lower-order

solution equations, that include the effects of higher-order disturbance

moments. Thus, there is a double infinity of possible approximate solutions: the

coefficients of a solution of any order may include the effects of disturbance

moments up to any order.

In the above 4th-order solution, to consider non-Gaussian disturbance

distributions that are symmetric about zero, we would need only to replace the

2nd and 4th moments, vec(Σ) and vec(Σ ⊗ Σ), with the relevant alternatives. For

example, without reference to any particular known distribution, we could fatten

the tails of the distribution by scaling up vec(Σ ⊗ Σ) in (3.41), while keeping

vec(Σ) constant everywhere. The solution equations must be rederived for

distributions that are not symmetric about zero.

4. NUMERICAL EXAMPLE.

This section illustrates the 4th-order solution equations derived in

Section 3. The example represents the optimizing behavior of a representative

firm in an industry, which maximizes an LQ approximation of its expected present

value, subject to its production function, internal adjustment costs, and output

demand and input supply conditions (Chen and Zadrozny, 2000a). The problem and

its linear feedback solution are given in the standard LQ form by the 2×2

matrices F0, G0, R0, S0, Q0, and P0, specifically,
18

(4.1) F0 = 





















.0000 .4571

.4056 .0000

, G0 = 





















.0000 .6935

.6587 .0000

, R0 = 





















−

−

44.33 44.43

44.43 44.33

,

 S0 = 





















− −

− −

.4351 .9350

.4349 .9338

, Q0 = 





















.8689 2.865

1.064 .8688

, P0 = 





















.1002 .1857

.1008 .1961

,

which imply

(4.2) Q = 





















.5189 2.164

.8912 .5189

, Φ = 





















.0695 .5858

.4720 .1292

, V = 





















1.037 3.566

1.256 1.037

.

To compute the 4th-order approximate solution of the LQEG problem, we also

assume that vec(Σ) = (1., .5, .5, 1.)T

 and evaluate the solution equations at x0

= ( .5 , .5 )T ≅ (.7071, .7071)T. Thus, we obtain the following values for the

coefficients of 4th-order approximate solution,

(4.3) g0 = .0043, ∇g1 = [.0021, .0059], ∇2

g0 = 





















.0081 .0241

.0049 .0119

,

 ∇3

g0

T

 = 



















−

.0608 .0608 .0784 .1213

.2663 .0356 .0435 .0610

,

 ∇4

g0

T

 = 





















.0622 .1185 .1179 .2334 .1149 .2233 .2244 .4988

.0390 .0637 .0642 .1110 .0699 .1196 .1217 .2341

.

These results give some sense of the sizes of the coefficients of the 4th-order

approximate solution.

5. CONCLUSION.
19

The paper has derived and illustrated recursions for computing any

multivariate Gaussian moment and equations for computing the 4th-order Taylorseries approximation of the objective function of the infinite-horizon LQEG

problem. There are several possible extensions of the paper. We could optimize

g(x) numerically with respect to the linear feedback matrix, P0. Of course, we

could also do this following Jacobson (1973), by solving an algebraic Riccati

equation. It would be interesting to derive and illustrate corresponding results

for a non-Gaussian disturbance distribution, say, the fatter-tailed multivariate

t distribution. Of course, for a non-Gaussian distribution, the optimal decision

rule will usually be nonlinear. In such a case, we should derive and solve

equations for simultaneously computing gˆ(x) and the approximate optimal

nonlinear decision rule, Pˆ (x), although this possibility would be significantly

more complicated than optimizing gˆ(x) numerically with respect to P0. In any of

these possibilities, we could consider the discounted problem.

The present perturbation method is different from that of Collard and

Juillard (2000) and Sims (2000), in which the stochastic problem is solved by

first perturbing the disturbance vector in the nonstochastic version of a

problem and, then, accounting for stochastic variation of the disturbances by

taking expectations of the nonstochastic perturbation equations. In our

experience, this method produces identical 1st- and 2nd-order approximate

solutions of stochastic problems, with equivalent effort, but, for higher-order

solutions of multivariate stochastic problems, produces intractable solution

equations. For this reason, we follow Fleming (1971) and Fleming and Souganidis

(1986) and account for stochastic variation in terms of the perturbation

parameter θ. Fleming’s method is also applied by Anderson and Hansen (1996),

Judd (1998), and Chen and Zadrozny (2000b).
20

APPENDIX A: DEFINITIONS AND RULES OF MATRIX DIFFERENTIATION.

A.1. DEFINITIONS OF MATRIX DERIVATIVES.

Let A(x) ∈ Dk

: Rn → Rp×q be a K-times differentiable p×q matrix function of

the n-vector x. A(x) could be a function of the matrix X ∈ Rk×m

, such that x =

vec(X), where vec(⋅) is the columnwise vectorization of a matrix. We consider

derivatives of elements of A with respect to elements of x in three forms: the ∂

or partial derivative form, the d or differential form, and the ∇ or gradient

form.

For k = 1, ..., K and i1, ..., ik ∈ {1, ..., n}, we define k

i1Lik ∂ A ∈ Rp×q

 by

(A.1) k

i1Lik ∂ A = 









































∂ ∂

∂

∂ ∂

∂

∂ ∂

∂

∂ ∂

∂

i1 ik

pq

k

ii ik

p1

k

i1 ik

1q

k

i1 ik

11

k

x x

A

x x

A

x x

A

x x

A

L L L

M M

L L L

,

as the partial derivative form of k-th order partial derivatives of the elements

of A with respect to i1 x , ..., ik x .

The differential form associated with (A.1) is

(A.2) dkA = ∑ = ∑ = ∂ n

ik 1 i1 ik

n

i1 1 L L A⋅d i1 x ⋅⋅⋅ d ik x ,

where the dxi’s are small (strictly, infinitesimal) increments to the elements

of x = (x1, ..., xn)T

.

The gradient form associated with (A.1) and (A.2) can now be built up

recursively, starting with k = 1. Because there is no generally accepted

terminology for higher-order derivatives of matrix-valued functions, such as

gradient and Hessian for scalar-valued functions or Jacobian for vector-valued

functions, we call the matrix representations of kth-order derivatives of matrix

functions "kth-order gradients."

Three rules are needed for vectorizing and permuting matrix elements.

First,
21

(A.3) vec(ABC) = [CT ⊗ A]vec(B),

where A, B, and C are matrices conformable to the matrix product ABC and ⊗

denotes the Kronecker matrix product (Magnus and Neudecker, 1988, p.30).

In sections 2 and 3, we apply (A.3) directly and inversely. "Directly"

means applying (A.3) from left to right, i.e., expressing the vectorization of

ABC as [CT ⊗ A]vec(B). "Inversely" means applying (A.3) from right to left,

i.e., expressing the unvectorization of [CT ⊗ A]vec(B) as a matrix whose

vectorization equals vec(ABC). We denote the inverse or "un" vectorization by

mat(⋅) and state its particular dimensions alongside.

Second,

(A.4) vec(A ⊗ B) = [In ⊗ (Kq,m ⊗ Ip)(Im ⊗ vec(B))]vec(A)

 = [(In ⊗ Kq,m)(vec(A) ⊗ Iq) ⊗ Ip]vec(B),

where A and B are m×n and p×q matrices, Ij denotes the j×j identity matrix, and

Kq,m denotes the qm×qm permutation matrix that maps the vectorization of a q×m

matrix to the vectorization of its transpose (Magnus and Neudecker, 1988, p.

48). The definition of Kq,m implies Kq,mT

 = (Kq,m)-1 = Km,q.

Third,

(A.5) (A ⊗ bT

)Kn,p = bT ⊗ A,

where A is m×n and bT

 is 1×p (Magnus and Neudecker, 1988, p. 47).

In B = Km,nA, where A and B are mn×p, Km,nA can be viewed as a permutation

operation on the elements of matrix A or as the matrix Km,n that premultiplies

the matrix A in the ordinary matrix product. In the first case, consider B =

Km,nA as bj = Km,naj, for j = 1, ..., p, where aj and bj are mn×1 columns of A and

B. For each j = 1, ..., p, consider aj as an m×n matrix, transpose this matrix,

and assign its column vectorization to bj. Thus, we construct B = Km,nA. In the

second case, consider Km,n = Km,nImn, where Imn is the mn×mn identity matrix. Km,n

on the right side of the equality is viewed as the permutation operator, as in

the first case, and Km,n on the left side is the result of permuting Imn. Thus,

we construct Km,n as a permutation of Imn.
22

We emphasize that throughout the paper a "vector" or the result of a

"vectorization" is a column vector.

A.2. REPRESENTATIONS OF MATRIX DERIVATIVES.

For k = 1, (A.1) and (A.2) become

(A.6) ∂kA = 









































∂

∂

∂

∂

∂

∂

∂

∂

i

pq

i

p1

i

1q

i

11

x

A

x

A

x

A

x

A

L

M M

L

 ,

(A.7) dA = ∑ =

n

i 1 ∂iA⋅dxi.

Note that vectorization, summation, and differentiation operations are

commutative, i.e., can be applied in any order. Therefore, we vectorize (A.7),

to obtain

(A.8) vec(dA) = [∂1vec(A), ..., ∂nvec(A)]dx,

where dx = (dx1, ..., dxn)T

, so that

(A.9) vec(dA) = ∇A⋅dx,

(A.10) ∇A = [∂1vec(A), ..., ∂nvec(A)].

Equations (A.9) and (A.10) relate the ∂, d, and ∇ forms of first-order

derivatives of A to each other.

To obtain analogues of (A.9) and (A.10) for k = 2, we differentiate them

to obtain

(A.11) vec(d2

A) = d(∇A)dx,

(A.12) d(∇A) = [d(vec(∂1A)), ..., d(vec(∂nA))]
23

 = ∑ =

n

j 1 [∂j(vec(∂1A)), ..., ∂j(vec(∂nA))]dxj.

Then, we vectorize (A.12) to obtain

(A.13) vec(d(∇A)) = ∑ =

n

j 1

























∂ ∂

∂ ∂

(vec( A))

(vec( A))

j n

j 1

M dxj = 





























∂ ∂ ∂ ∂

∂ ∂ ∂ ∂

(vec( A)) (vec( A))

(vec( A)) (vec( A))

1 n n n

1 1 n 1

L

M M

L

dx.

Then, because vec(∇A) = 

























∂

∂

(vec(A))

(vec(A))

n

1

M = 

























∂

∂

vec( A)

vec( A)

n

1

M , we obtain

(A.14) vec(d(∇A)) = [∂1vec(∇A), ..., ∂nvec(∇A)]dx.

Continuing in this manner for k = 2, ..., K, we obtain

(A.15) vec(d(∇k-1A)) = ∇k

A⋅dx,

(A.16) vec(dk

A) = [(Πk-1⊗dxT) ⊗ Ipq]∇kA⋅dx,

where Πk-1⊗dxT

 denotes k-1 successive Kronecker products of dxT, and

(A.17) ∇k

A = [∂1vec(∇k-1A), ..., ∂nvec(∇k-1A)].

Applied for k = 1, ..., K, (A.17) recusively organizes gradient form derivatives

of A up to order K as matrices. Basically, ∇kA is the Jacobian matrix of the

vectorization of ∇k-1A.

In this gradient representation of matrix derivatives, the K-term Taylorseries approximation of A(x) at x = x0 is

(A.18) vec( Aˆ (x)) = ∑ =

K

k 0 (1/k!)vec(∇k

A0)T{[Πk⊗(x-x0)] ⊗ Ipq},

such that ∇k

A0 = ∇kA(x0), for k ≥ 1, and ∇0A0 = A(x0).
24

A.3. DIFFERENTIATION RULES.

Let A(x) ∈ D: Rn → Rp

 and B(y) ∈ D: Rp → Rq be differentiable vector

functions (or vectorizations of matrix functions). Let C(x) ∈ D: Rn → Rq

 be the

differentiable composite vector function C(x) = B(A(x)). Then, the gradient form

of derivatives of C(x) is given by the chain rule of differentiation,

(A.19) ∇C(x) = ∇B(A)⋅∇A(x).

Let A(x) ∈ D: Rn → Rp×q

 and B(x) ∈ D: Rn → Rp×q be differentiable matrix

functions conformable to the ordinary matrix product C(x) = A(x)⋅B(x). Then, the

differential form of derivatives of C(x) is given by the product rule of

differentiation,

(A.20) dC(x) = dA(x)⋅B(x) + A(x)⋅dB(x).

Rules (A.19) and (A.20) are quickly proved by elementwise application of

the scalar chain rule of differentiation and the scalar product rule of

differentiation.
25

APPENDIX B: INPUT DERIVATIVES FOR PERTURBATION SOLUTION.

This appendix derives equations for computing the first to fourth gradient

matrices ∇f(x), ..., ∇4

f(x) of f(x) = exp[-(1/2)xTQx], which are the inputs for

computing the 4th-order approximation of the stochastic solution function, g(x).

As in Section 3, without further reference, we vectorize scalars, vectors, and

matrices using rule (A.3), differentiate expressions using the chain and product

rules of differentiation, (A.19) and (A.20), and repeatedly use the fact that

f(x) is scalar valued. However, we do reference vectorization rule (A.4) when we

use it. In the example in section 4, f(x), ..., ∇4

f(x) are evaluated at x0.

Differentiating f(x), we obtain

(B.1) df(x) = -f(x)xT

Qdx.

Then, noting that df(x) = ∇f(x)dx and dropping dx, we obtain

(B.2) ∇f(x) = -f(x)xT

Q.

Differentiating (B.2), we obtain

(B.3) d(∇f(x)) = -df(x)xT

Q - f(x)dxTQ.

Vectorizing (B.3) and eliminating df(x) using (B.1), we obtain vec[d(∇f(x))] =

f(x)(QxxT

Q + Q)dx. Then, noting that vec[d(∇f(x))] = ∇2f(x)dx and dropping dx,

we obtain

(B.4) ∇2

f(x) = f(x)(-Q + QxxTQ).

Differentiating (B.4), we obtain

(B.5) d(∇2

f(x)) = df(x)(-Q + QxxTQ) + f(x)Q(dxxT + xdxT)Q.

Vectorizing (B.5), we obtain vec[d(∇2

f(x))] = vec(-Q + QxxTQ)df(x) + f(x)(Qx ⊗

Q)dx + f(x)(Q ⊗ Qx)dx. Then, eliminating df(x) using (B.1), noting that

vec[d(∇2

f(x))] = ∇3f(x)dx, and dropping dx, we obtain

(B.6) ∇3

f(x) = f(x)[vec(Q - QxxTQ)xTQ + (Qx ⊗ Q) + (Q ⊗ Qx)].
26

Differentiating (B.6), we obtain

(B.7) d(∇3

f(x)) = df(x)[vec(Q - QxxTQ)xTQ + (Qx ⊗ Q) + (Q ⊗ Qx)]

 + f(x){-vec[Q(dxxT

 + xdxT)Q)]xTQ + vec(Q - QxxTQ)dxTQ

 + (Qdx ⊗ Q) + (Q ⊗ Qdx)}.

Vectorizing (B.7), we obtain

(B.8) vec[d(∇3

f(x))] = vec[vec(Q - QxxTQ)xTQ + (Qx ⊗ Q) + (Q ⊗ Qx)]df(x)

 - f(x){(Qx ⊗ 2 nI )[(Qx ⊗ Q) + (Q ⊗ Qx)]dx

 + [Q ⊗ vec(Q - QxxT

Q)]dx

 + vec(Qdx ⊗ Q) + vec(Q ⊗ Qdx)},

where 2 nI denotes the n2

×n2 identity matrix. Vectorization rule (A.4) implies

(B.9) vec(Qdx ⊗ Q) = [Kn,n ⊗ In][In ⊗ vec(Q)]Qdx,

(B.10) vec(Q ⊗ Qdx) = [vec(Q) ⊗ In]Qdx,

because I1 = scalar 1 and Kn1 = K1n = In. Then, using (B.1), (B.9), and (B.10),

noting that vec[d(∇3

f(x))] = ∇4f(x)dx, and dropping dx, (B.8) becomes

(B.11) ∇4

f(x) = f(x){vec[vec(Q - QxxTQ)xTQ + (Qx ⊗ Q) + (Q ⊗ Qx)]xTQ

 - (Qx ⊗ 2 nI )[(Qx ⊗ Q) + (Q ⊗ Qx)] + [Q ⊗ vec(Q - QxxT

Q)]

 + [Kn,n ⊗ In][In ⊗ vec(Q)]Q + [vec(Q) ⊗ In]Q}.

REFERENCES.
27

Anderson, E.W. and L.P. Hansen (1996), "Perturbation Methods for Risk-Sensitive

Economies," working paper, Department of Economics, University of Chicago,

Chicago, IL.

Chen, B. and P.A. Zadrozny (2000a), "Estimated U.S. Manufacturing Capital and

Productivity Based on an Estimated Dynamic Economic Model," working paper,

Department of Economics, Rutgers University, Camden, NJ.

Chen, B. and P.A. Zadrozny (2000b), "Perturbation Solution of Nonlinear Rational

Expectations Models," working paper, presented at the Computational Economics

segment of the SITE 2000 summer workshop, Department of Economics, Stanford

University, Stanford, CA.

Collard, F. and M. Juillard (2000), "Accuracy of Stochastic Perturbation

Methods: The Case of Asset Pricing Models," Journal of Economic Dynamics and

Control, forthcoming.

Fleming, W.H. (1971), "Stochastic Control for Small Noise Intensities," SIAM

Journal of Control 9: 473-517.

Fleming, W.H. and P.E. Souganidis (1986), "Asymptotic Series and the Method of

Vanishing Viscosity," Indiana University Mathematics Journal 35: 425-447.

Hansen, L.P. and T.J. Sargent (1995), "Discounted Linear Exponential Quadratic

Gaussian Control," IEEE Transactions on Automatic Control AC-40: 968-971.

Jacobson, D.H. (1973), "Optimal Stochastic Linear Systems with Exponential

Performance Criteria and their Relation to Deterministic Differential Games,"

IEEE Transactions on Automatic Control AC-18: 124-131.

Judd, K.L. (1998), Numerical Methods in Economics, Cambridge, MA: MIT Press.

Karp, L.S. (1985), "Higher Moments in the Linear-Quadratic-Gaussian Problem,"

Journal of Economic Dynamics and Control 9: 41-54.

Laub, A.J. (1979), "A Schur Method for Solving Algebraic Riccati Equations," IEEE

Transactions on Automatic Control 24: 913-921.

MacRae, E.C. (1974), "Matrix Derivatives with an Application to an Adaptive

Linear Decision Problem," Annals of Statistics 2: 337-346.

Magnus, J.R. and H. Neudecker (1988), Matrix Differential Calculus with

Applications in Statistics and Econometrics, New York, NY: J. Wiley & Sons.

Sims, C.A. (2000), "Second Order Solution of Discrete Time Dynamic Equilibrium

Models," working paper, Department of Economics, Princeton University,

Princeton, NJ.