Mathematica Notebook for "Asymptotic Methods for Asset Market Equilibrium Analysis" by Judd and Guu

The following Mathematica program computes asset market equilibrium for two investors, one safe asset, and one or two risky assets.

The first command turns off annoying spelling queries.

Off[General::spell1]

The initial setup

Economic Model

Define the two risky assets' returns, Z and Y, in terms of zero-mean random variables z and y. ϵ is the scaling parameter, equal to the standard deviation in the ϵ-economy. Z is an asset with mean return R and variance (ϵ^2), and is in positive net supply. It represents equity. To keep the formulas simple, we assume that the return on bonds, R, is 1. This is equivalent to assuming that the bond in the second period is the numeraire.
	These facts imply that Z can be decomposed in the following manner:

Z = 1 + ϵ z;

Y is the derivative asset and has zero net supply. It has mean return ϵ μY, is partially correlated with Z through the ϵ α z term, and has an orthogonal component ϵ y. These facts imply that Y can be represented as

Y = ϵ μY   +  ϵ α z  + ϵ  y;

The price of Z (Y) is p (q). We parametrize them in terms of the scaling parameter ϵ, the mean of Y,  μY, and the premia of Z and Y, denoted Π and Ψ:

p  =  1  - (ϵ^2)  Π ;  
    q  = ϵ  μY  -  (ϵ^2) Ψ ;

Note: The notation in the notebook corresponds to the notation in the paper except for a minor change. In Mathematica, the letter π is reserved for 3.14159... and so we could not use it as the risk premium for Z. Therefore, we use Π instead. To maintain the symmetry we also let Ψ denote the premium for Y.

(Subscript[W, i]) is final wealth of type i investor. θ is type 1 demand for equity. θ(Subscript[e, i]) is type i endowment of equity and (Subscript[Be, i]) is type i endowment of bonds.  ϕ is type 1's demand for Y; - ϕ is type 2 demand. In equilibrium, type 2 investors will hold θ(Subscript[e, 1])+θ(Subscript[e, 2])-θ shares of stock. ( Mathematica prevented us from using (Subscript[θ, i]^e) ((Subscript[B, i]^e)) to denote type i endowment of the risky asset (bond).)

θ(Subscript[e, 1]) = θ(Subscript[e, 1]); θ(Subscript[e, 2]) = θ(Subscript[e, 2]);

(Subscript[W, 1])=((Subscript[Be, 1])+ p θ(Subscript[e, 1])-p θ- q ϕ)+θ*Z+ϕ Y
(Subscript[W, 2])=((Subscript[Be, 2])+ p θ(Subscript[e, 2])-p (θ(Subscript[e, 1]) + θ(Subscript[e, 2]) - θ)- q (-ϕ))+(θ(Subscript[e, 1]) + θ(Subscript[e, 2]) - θ)*Z+(-ϕ) Y;
(Subscript[W, 1]) = Expand[(Subscript[W, 1])];
(Subscript[W, 2]) = 
        Expand[(Subscript[W, 2])];

We compute the four first-order conditions for the two investors and the two risky assets. The bond demand is determined by substracting risky asset demand from initial wealth. We divide each first-order condition by ϵ to eliminate one degree of degeneracy. FOCij is the first-order condition of type i investor with respect to risky asset j, where j = 1 is equity and j = 2 refers to the derivative.

FOC11 = (1/ϵ) D[(Subscript[u, 1])[(Subscript[W, 1])], θ] // 
        Simplify
FOC12 = (1/ϵ) D[(Subscript[u, 1])[(Subscript[W, 1])], ϕ] // Simplify
FOC21 = (1/ϵ) D[(Subscript[u, 2])[(Subscript[W, 2])], θ] // 
        Simplify
FOC22 = (1/ϵ) D[(Subscript[u, 2])[(Subscript[W, 2])], ϕ] // 
        Simplify

Even though we have divided by ϵ, these first-order conditions describe the situation at ϵ=0.
These four first-order conditions define equilibrium. Define G to be the vector of first-order conditions.

Moment substitutions and other substitutions

We next define lists of substitutions that allow us to compute the moments of z and y. The doubly subscripted term (Subscript[μ, m, n]) denotes the expectation of (z^m) (y^n). We often need to compute expectations of our Taylor series expansions. We could invoke Mathematica's Integrate command, but we can do better by searching for specific terms in our polynomials and replacing them with their integrals. For example, whenever we see an isolated z in an expression we know that its integral just replaces the z term by the mean of z, which is denoted as (Subscript[μ, 1, 0]). This works better than integration in Mathematica since we know that our function is a polynomial whereas Mathematica needs to spend time to ascertain that fact. The following replacement rules implement this approach to integration.

ownmoments  =  {z -> (Subscript[μ, 1, 0]), 
        (z^m)_ ->  (Subscript[μ, m, 0]), y -> (Subscript[μ, 0, 1]), 
        (y^m)_ ->  (Subscript[μ, 0, m])}
crossmoments  =  {(z^m)_ (y^n)_ ->  (Subscript[μ, m, n]), 
        z (y^n)_ ->  (Subscript[μ, 1, n]), (z^m)_ 
            y ->  (Subscript[μ, m, 1])}

Our definition of Y decomposed it into its mean, it covariance with Z, and the orthogonal residual. Without loss of generality, we set the mean of z and y to zero.

zeromean  =  {(Subscript[μ, 0, 1]) -> 0, (Subscript[μ, 1, 0]) -> 
            0};

We make z and y orthogonal and set the variance of z equal to 1.

orthogonal  =  {(Subscript[μ, 1, 1]) -> 0, (Subscript[μ, 2, 0]) -> 
            1};

We now gather all the substitions into one list.

moments  =  
      Join[ownmoments, crossmoments, zeromean,  orthogonal]

The following substitutions denote the case where there is no derivative security. They say that all cross-moments between z and y are zero. Note that y must be nonzero since we set its variance to 1. However, if all of its moments are uncorrelated with all moments of the endowment then it is a pure gamble and will not be traded in equilibrium.

NoDeriv = 
      Table[(Subscript[μ, m, n]) -> 0, {m, 1, 6}, {n, 1, 6}] // 
        Flatten

We define a convenient list of simplifying substitutions which Mathematica did not automatically execute.

sub1={(1/2)(2 (Subscript[Be, 1])+2 θ(Subscript[e, 1]))-> ((Subscript[Be, 1])+θ(Subscript[e, 1])),(1/2)(2 (Subscript[Be, 2])+2 θ(Subscript[e, 2]))->((Subscript[Be, 2])+θ(Subscript[e, 2]))}

Solvability Condition in Theorem 4

We check for the nonsingularity of the solvability matrix, (Subscript[G, ϵ, Λ]) where Λ=(θ,ϕ,Π,Ψ). First define the solvability matrix:

(Subscript[G, ϵ, Λ])=
{{D[FOC11, {ϵ, 1}, θ], 
            D[FOC11, {ϵ, 1}, ϕ], 
            D[FOC11, {ϵ, 1}, Ψ], 
            D[FOC11, {ϵ, 1}, Π]}, 
{D[
              FOC12, {ϵ, 1}, θ], 
            D[FOC12, {ϵ, 1}, ϕ], 
            D[FOC12, {ϵ, 1}, Ψ], 
            D[FOC12, {ϵ, 1}, Π]}, 
{D[
              FOC21, {ϵ, 1}, θ], 
            D[FOC21, {ϵ, 1}, ϕ], 
            D[FOC21, {ϵ, 1}, Ψ], 
            D[FOC21, {ϵ, 1}, Π]}, 
{D[
              FOC22, {ϵ, 1}, θ], 
            D[FOC22, {ϵ, 1}, ϕ], 
            D[FOC22, {ϵ, 1}, Ψ], 
            D[FOC22, {ϵ, 1}, Π]}};

Evaluate (Subscript[G, ϵ, Λ]) at ϵ=0 and display it in matrix form:

(Subscript[G, ϵ, Λ]) = 
        (Subscript[G, ϵ, Λ]) /. {ϵ -> 0};
(Subscript[G, ϵ, Λ]) // MatrixForm

Integrate the elements of (Subscript[G, ϵ, Λ]) by replacing powers of z and y with their moments.

ExpG  =  
        Expand[(Subscript[G, ϵ, Λ])] //. 
          moments;
ExpG // MatrixForm

Compute the determinant of E {(Subscript[G, ϵ, Λ])}.

(Det[ExpG] // Simplify) //. moments

As long as (Subscript[μ, 0, 2]) is not zero, this determinant is strictly positive and (Subscript[G, ϵ, Λ]) is nonsingular

Bifurcation equations: series construction

We construct substitutions that will define the equilibrium mappings from the scaling parameter ϵ to the equilibrium values of (θ,ϕ,Π,Ψ).

PortfolioSubs  =  {θ -> θ[ϵ], ϕ -> ϕ[ϵ]}

PremiaSubs  =  {Π -> Π[ϵ], Ψ -> Ψ[ϵ]}

To construct the equilibrium equations we take the first-order conditions and replace the the portfolio variables, θ and ϕ, and the risk premium variables, Π and Ψ,  with their equilibrium maps.

EQM11  =  FOC11 /. PortfolioSubs /. PremiaSubs
EQM21 =  FOC21 /. PortfolioSubs /. PremiaSubs
EQM12  =  FOC12 /. PortfolioSubs /. PremiaSubs
EQM22 =  FOC22 /. PortfolioSubs /. PremiaSubs

The equilibrium equations are really the expectations 0=E{EQMij}. We will execute the integrations later, after the Taylor series expansions, since it is permissible to interchange the integration and differentiation operations.

We define abbrevations for the derivatives; this helps make the formulas more understandable than the ones automatically generated by Mathematica.

UtilDerivs = {Derivative[n_][(Subscript[u, i])_][(Subscript[Be, i])_ + θ(Subscript[e, i])_] -> 
          (Subscript[A, i])[n]}

We want to replace derivatives of the utility function, the (Subscript[A, i])[m] terms, with risk tolerance, (Subscript[τ, i]), skew tolerance, (Subscript[ρ, i]), and kurtosis tolerance, (Subscript[κ, i]). The next substitution rule does that.

UtilParams  =  {(Subscript[A, i])_[2] -> -(Subscript[A, i])[1]/(Subscript[τ, i]), 
        (Subscript[A, i])_[3] -> ((2 (Subscript[ρ, i]) (Subscript[A, i])[1])/((Subscript[τ, i])^2)), 
        (Subscript[A, i])_[4] -> ((3 (Subscript[κ, i]) (Subscript[A, i])[1])/((Subscript[τ, i])^3)), 
        (Subscript[A, i])_[1] ->  1}

It is often useful to refer to the social risk tolerance, T, and express type i risk tolerance as a share, (Subscript[ν, i]), of total risk tolerance. The following replacement rule allows us to do that.

taureps  =  {(Subscript[τ, 1]) + (Subscript[τ, 2]) -> 
            T, (Subscript[τ, i])_ -> (Subscript[υ, i]) 
              T, (Subscript[υ, 1]) + (Subscript[υ, 2]) -> 1};

We now compute the power series expansions of the four equilibrium conditions. We fully expand each series so that the various powers of z and y are gathered together. This is necessary for our integration approach to work.

EqmPow11  =  (Normal[Series[EQM11 , {ϵ, 0, 6}]] /. 
              sub1) //. UtilDerivs;
EqmPow11  =  Expand[EqmPow11];
EqmPow12  =  (Normal[Series[EQM12 , {ϵ, 0, 6}]] /. 
              sub1) //. UtilDerivs;
EqmPow12  =  Expand[EqmPow12];
EqmPow21  =  (Normal[Series[EQM21 , {ϵ, 0, 6}]] /. 
              sub1) //. UtilDerivs;
EqmPow21  =  Expand[EqmPow21];
EqmPow22  =  (Normal[Series[EQM22 , {ϵ, 0, 6}]] /. 
              sub1) //. UtilDerivs;
EqmPow22  =  Expand[EqmPow22];

We now take the expectation of each equilibrium equation's power series. Since z and y are the only random variables, we replace each power of z and y and each crossproduct by the appropriate moment.

EqmPow11  =  EqmPow11 //. moments;
EqmPow12  =  EqmPow12 //. moments;
EqmPow21  =  EqmPow21 //. moments;
EqmPow22  =  EqmPow22 //. moments;

The final step in constructing the equilibrium expressions is to collect terms of like powers of ϵ and list the coefficients of each power of ϵ. This puts the equations in the proper arrangement for solving the problem.

Eqns11  =  CoefficientList[EqmPow11, ϵ];
Eqns12  =  CoefficientList[EqmPow12, ϵ];
Eqns21  =  CoefficientList[EqmPow21, ϵ];
Eqns22  =  CoefficientList[EqmPow22, ϵ];

Solve the individual equations in sequence

Equation 1:

The coefficients of the (ϵ^0) components of the equilibrium power series should be zero. We now check what they are.

{Eqns11[[1]],  Eqns21[[1]],  Eqns12[[1]],  
      Eqns22[[1]]}

This shows that we can continue. If this were not a vector of zeroes then we would know that our parameterization did not fulfill the necessary conditions of the Bifurcation Theorem.

Equation 2:

We compute the bifurcation point by choosing (θ[0],ϕ[0],Ψ[0],Π[0]) so that the ϵ components of the equilibrium equation expansion are zero. We first list the equations

Eq1 = Eqns11[[2]] // Simplify // Expand
Eq2 = Eqns21[[2]] // Simplify // Expand
Eq3 = Eqns12[[2]] // Simplify // Expand
Eq4 = Eqns22[[2]] // Simplify // Expand

We next solve for bifurcation point.

BifPt =  Solve[{Eq1 == 0,  Eq2 == 0,  
            Eq3 == 0,  
            Eq4 == 0}, { θ[0], ϕ[0], Ψ[
              0], Π[0]}][[1]]

We simplify the expression, bringing together common factors.

BifPt  = BifPt // Simplify

We next simplify by using the equilibrium conditions. Essentially, we want to replace all occurrences of θ(Subscript[e, 1]) + θ(Subscript[e, 2]) terms, the aggregate endowment of the risky asset, with its value, 1. In geneal, we let Θ be the total endowment of the risky asset; we set Θ=1. EqmSubs defines various substitutions that express the identity θ(Subscript[e, 2]) + θ(Subscript[e, 1]) = Θ, andvec = -  allow us to simplify various expressions.

Θ = 1;  
    EqmSubs  =  {θ(Subscript[e, 2]) + θ(Subscript[e, 1]) -> Θ, θ(Subscript[e, 2]) + θ(Subscript[e, 1]) - Θ -> 
          0, Θ - θ(Subscript[e, 2]) - θ(Subscript[e, 1]) -> 
          0}

Simplify BifPt by uskng the substitutions in EqmSubs

BifPt  = BifPt //. EqmSubs

We do not like the utility derivatives,  (Subscript[A, i])[j]. We apply substitutions contained in UtilParams, defined above, that replace utility derivatives with indices such as risk tolerance.

BifPt  =  BifPt //. UtilParams

We need to simplify this expression for BifPt. The Simplify command can handle this for us. (In general, one must be careful using the Simplify command since it can often take a long time to find the desired simplification or it will find a simplification other than the one you want.)

BifPt  =  BifPt // Simplify

We now have an intuitive expression for the bifurcation point. The procedure was direct: solve for the bifurcation point in terms of utility derivatives and endowments, and then apply various simplifications to transform the solution into an expression involving elasticities and shares. Below we will apply similar sequences of simplifications without explanation.

The bifurcation point conforms to basic intutions. The risk premium, Π[0], equals the inverse of the total risk tolerance,  (1/((Subscript[τ, 1]) + (Subscript[τ, 2]))). Type 1 holdings of equity equals it share of social risk tolerance. The risk premium of the derivative asset Y, Ψ[0], equals (α/((Subscript[τ, 1]) + (Subscript[τ, 2])))which is the product of the risk premium of equity and the covariance of Y with z. (This is equation (32) in the paper). The most surprising result is that the asymptotic holding of the derivative asset is zero. This does not mean that there is no demand for the derivative; it means that demand is smaller than order 0. We will see below that demand is of order ϵ.

Equation 3:

We next compute (θ'[0],ϕ'[0],Ψ'[0],Π'[0]) by setting the (ϵ^2) components of the equilibrium equations Taylor series expansions equal to zero. We immediately make the substitutions stated in UtilParams.

Eq1  =  Eqns11[[3]] //. UtilParams // Expand;
Eq2  =  Eqns21[[3]] //. UtilParams // Expand;
Eq3  =  Eqns12[[3]] //. UtilParams // Expand;
Eq4  =  Eqns22[[3]] //. UtilParams // Expand;

Here is the first equation; the rest are similar.

Eq1

We now solve for the unknowns, θ'[0], ϕ'[0], Ψ'[0], and Π'[0]. We suppress printing the first set of results since they take up too much space. The substitutions we make below will produce comprehensible expressions.

Sol2  =  Solve[{Eq1 == 0,  Eq2 == 0,  
              Eq3 == 0,  
              Eq4 == 0}, { θ'[0], ϕ'[
                0], Ψ'[0], Π'[
                0]}][[1]];

The solution involves , θ[0], ϕ[0], Ψ[0], and Π[0]. We now replace these terms with the values found in our bifurcation point solution, BifPt.

Sol2  =  Sol2 /. BifPt;

Note: we could have made these substitutions when we defined Eq1, Eq2, etc, but that would not be wise. If we had done that then the expressions in Eq1, Eq2, etc., would have been larger, making more work for Mathematica. This would not have been a major problem for this set of equations, but becomes important when we move to higher-order terms where the solutions become more complex.

We next simplify the result which produces a fairly compact form.

Sol2  =  Sol2 // Simplify;
Sol2 // TableForm

We can further simplify it by applying the equilibrium substitutions.

Sol2  =  Sol2 /. EqmSubs // Simplify;
Sol2 // TableForm

The solution for ϕ'[0] in Sol2 proves the assertions in Theorem 9.

If there is no derivative security, then the solution is

Sol2NoY = Sol2 //. NoDeriv

This proves Theorem 8. The solution for θ'[0] in Sol2NoY proves the assertion in equation (26) and the solution for Π'[0]in Sol2NoY proves the assertion in equation (27).

Compute the change in equity demand from introduction of Y

Notice that the presence of Y has no impact on the equity premium derivative, Π'[0], but it does affect equity demand, θ'[0]. So, we now determine that effect.

If there are two assets, then θ'[0] is

TwoAsset  =  θ'[0] /. BifPt /. Sol2

If α=0 then the derivative asset has no impact on equilibrium and it is as if it did not exist. The substitutions in NoDeriv (defined above) state that the derivative is uncorrelated with equity at all moments, representing the case where there is no derivative. Therefore, equity demand in its absence is

SingleAsset  =  TwoAsset /. NoDeriv

The change in equity demand is the difference:

TwoAsset  -  SingleAsset // Simplify

This is the result reported in Theorem 10.

Equation 4:

We next compute (θ''[0],ϕ''[0],Ψ''[0],Π''[0]) by setting the (ϵ^3) components of the equilibrium equations equal to zero.

Eq1  =  Eqns11[[4]] //. UtilParams /. Sol2 /. 
              BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])} // Expand;
Eq2  =  Eqns21[[4]] //. UtilParams /. Sol2 /. 
              BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])} // Expand;
Eq3  =  Eqns12[[4]] //. UtilParams /. Sol2 /. 
              BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])} // Expand;
Eq4  =  Eqns22[[4]] //. UtilParams /. Sol2 /. 
              BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])} // Expand;

Here is the first equation; the rest are similar.

Eq1

Solve for the second derivatives

Sol3  =  Solve[{Eq1 == 0,  Eq2 == 0,  
              Eq3 == 0,  
              Eq4 == 0}, { θ''[0], ϕ''[
                0], Ψ''[0], Π''[
                0]}][[1]];

Sol3 = Sol3 // Simplify;
Sol3

Change in equity risk premium

The presence of Y does affect Π''[0]. pi2 expresses this term.

pi2  =  Π''[0] /. Sol3 // Together // 
        Expand

Define the price function

PriceZ[ϵ_]  =  
      1  - (ϵ^2)  Π[ϵ]

Compute the Taylor series for the price function

TaylorPriceZ =  
      Normal[Series[PriceZ[ϵ] , {ϵ, 0, 4}]]

TaylorPriceTwoAssetZ is the second-order expansion for the price of Z when there are two assets.

TaylorPriceTwoAssetZ = TaylorPriceZ //. UtilParams /. Sol3 /. 
              Sol2 /. 
            BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])} // 
        Simplify

TaylorPriceOneAssetZ is the second-order expansion for the price of Z when there is just Z. It equals TaylorPriceTwoAssetZ with the moments for Y and cross-moments between Z and Y zeroed out by the substitutions in NoDeriv.

TaylorPriceOneAssetZ = 
      TaylorPriceTwoAssetZ //. NoDeriv // Simplify

We now compute the difference to determine the impact of the derivative asset on the price of Z.

TaylorPriceTwoAssetZ - TaylorPriceOneAssetZ // Together // 
      Simplify

This is the result in Theorem 11.

Equation 5:

We next compute (θ'''[0],ϕ'''[0],Ψ'''[0],Π'''[0]) by setting the (ϵ^4) components equal to zero. We do this since these third derivatives may play a role in the utility analysis below.

Eq1  =  Eqns11[[5]] //. UtilParams /. Sol3 /. 
              Sol2 /. 
            BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])};
Eq2  =  Eqns21[[5]] //. UtilParams /. Sol3 /. 
              Sol2 /. 
            BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])};
Eq3  =  Eqns12[[5]] //. UtilParams /. Sol3 /. 
              Sol2 /. 
            BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])};
Eq4  =  Eqns22[[5]] //. UtilParams /. Sol3 /. 
              Sol2 /. 
            BifPt /. {θ(Subscript[e, 2]) -> Θ - θ(Subscript[e, 1])};

We display the first equation as an example

Eq1 // Simplify

We now solve for the third derivatives. We do not display any of the results since they are very long and difficult to interpret.

Sol4  =  Solve[{Eq1 == 0,  Eq2 == 0,  
              Eq3 == 0,  
              Eq4 == 0}, { θ'''[0], ϕ'''[
                0], Ψ'''[
                0], Π'''[0]}][[1]];

Utility Expansion

We want to evaluate the impact of the derivative security Y on utility. This is necessary to derive the results in Section 6.4.

Utility

(Subscript[W, 1]) is the final wealth and consumption of a type 1 agent. Express it in terms of the random variables z and y, and portfolio holdings.

(Subscript[W, 1]) = Expand[(Subscript[W, 1])]

Define equilibrium expected utility for type 1 agents by substituting the equilibrium functions for θ, ϕ, Π, and Ψ into our expression for final consumption. The result is equilibrium utility as a function of ϵ.

U1 = (Subscript[u, 1])[(Subscript[W, 1])] /. PortfolioSubs /. PremiaSubs

Compute the degree 5 Taylor series of the utility of type 1 agents in terms of ϵ and call it U1pow.

U1pow = 
        Normal[Series[U1 , {ϵ, 0, 5}]] //. UtilDerivs;

Expand it so that products of z and y are collected

U1pow  =  Expand[U1pow];

Compute expected utility by replacing instances of powers of z and y with their moment expressions

U1pow  =  U1pow //. moments;

Replace θ, ϕ, Π, and Ψ and their derivatives with the solutions derived above.

U1pow  =  U1pow /. Sol4 /. Sol3 /. Sol2 /. BifPt //. 
          UtilDerivs;

Express in terms of τ, ρ, and κ.

U1pow  =  U1pow //. UtilParams;

Combine like powers of ϵ and put the coefficients in the list U1ϵpows.

U1ϵpows  =  
        CoefficientList[U1pow, ϵ];

Display the coefficients of 1, ϵ, (ϵ^2), and (ϵ^3)

Table[U1ϵpows[[i]], {i, 1, 4}] // Simplify

This list shows that the moments of Y and its cross-moments with Z have no impact on utility up to the order (ϵ^3). Therefore, we move to the fourth-order terms to determine the utility impact of the derivative security.

Utility difference - order 4

Compute the utility contribution of the new asset. The utility difference is the utility with a nontrivial Y minus the utility with trivial Y, that is, a Y with zero co-moments with z.

Util4 is the degree four term in the utility expansion with two assets and equals

Util4  =  U1ϵpows[[5]] // Together

This term is very complex and difficult to interpret. However, since we are only interested in the impact of the new security on utility, we do not need to understand all of Util4. We only need to understand how the new security affects Util4. Examination of Util4 shows that (Subscript[μ, 2, 1]) is the only moment involving the new asset in the expression Util4. The degree four term in the utility expansion for the case with only one asset is computed by evaluating Util4 with (Subscript[μ, 2, 1]) set equal to zero.

(Util4 /. { (Subscript[μ, 2, 1]) -> 0})

We take the difference in the past two expressions to compute the impact of the new security.

UtilDiff = (Util4  - (Util4 /. { (Subscript[μ, 2, 1]) -> 0})) // 
        Expand

This is a much simpler expression. We simplify the difference to arrive at the result in the paper.

UtilDiff = Simplify[UtilDiff]

Use the identity -1 + θ(Subscript[e, 1]) - θ(Subscript[e, 2]):

UtilDiff = 
      UtilDiff /. {-1 + θ(Subscript[e, 1]) -> -θ(Subscript[e, 2])}

The expression in Theorem 12 is

Thm12 = ((((((Subscript[ρ, 1]) - (Subscript[ρ, 2])))^2) ((Subscript[τ, 1])^2) ((Subscript[τ, 2])^2)  )/(2 ((((Subscript[τ, 1]) + (Subscript[τ, 2])))^6) )) (4 (θ((Subscript[e, 1])/( (Subscript[τ, 1]))) - θ((Subscript[e, 2])/(Subscript[τ, 2]))) + 
            (1/( (Subscript[τ, 1])))) ((  (Subscript[μ, 2, 1]^2))/( (Subscript[μ, 0, 2])))

UtilDiff - Thm12 // Simplify

Therefore, our UtilDiff expression is equivalent to the one for utility change in Theorem 12.