Stoch Proc - Lab/Tutorial 5 Solutions

Scott Alexander

Q1

(a)

Theoretical variance reduction

f[x_]:=1/7; (* define density f *)
g[x_]:=21x (E^(-x/2)); (* define function g *)
q[x_]:=(Sqrt[x]); (* define control variate q *)
GG=(∫_1^8)g[x] f[x]\[DifferentialD]x; //Expand(* 𝔼[g(X)] *)
QQ=(∫_1^8)q[x] f[x]\[DifferentialD]x; (* 𝔼[q(x)] *)
varg=(∫_1^8)((g[x])^2) f[x]\[DifferentialD]x-(GG^2); (* var(g(X)) *)
varq=(∫_1^8)((q[x])^2) f[x]\[DifferentialD]x-(QQ^2); (* var(q(X)) *)
covgq=(∫_1^8)g[x]q[x] f[x]\[DifferentialD]x-GG QQ; (* cov(g(X),q(X)) *)
corgq=(covgq/(Sqrt[varg varq]));  (* corr(g(X),q(X)) *)
1-(corgq^2)//N (* ratio var((Subscript[(OverTilde[G]), n])) to var((Subscript[G, n])) *)

Using Mathematica

n=(10^5); (* sample size *)
g[x_]:=21x (E^(-x/2)); (* define function g *)
q[x_]:=(Sqrt[x]); (* define control variate q *)
QQ=Expectation[q[x],x\[Distributed]UniformDistribution[{1,8}]];   (* E[q(x)] *)

rv=RandomReal[UniformDistribution[{1,8}],n]; (* generate U(1,8) RVs *)
grv=g[rv]; (* apply g to RVs *)
qrv=q[rv]; (* apply q to RVs *)
aStar=Covariance[grv,qrv]/Variance[qrv]   (* calculate optimum a *)
rv=grv-aStar(qrv-QQ);    (* apply estimator to RVs *)

GnTilde=Mean[rv]; (* calculate mean of estimator *)
varGnTilde=Variance[rv]/n;  (* calculate sample variance of estimator *)
sdGnTilde=(Sqrt[varGnTilde]); (* calculate sample standard error *)
{GnTilde,varGnTilde,sdGnTilde}//N (* display results *)

Using R

See R code file.

(b)

Theoretical variance reduction

μX=Mean[UniformDistribution[{1,8}]];  (* mean of U(1,8) *)
g[x_]:=21x (E^(-x/2)); (* define function g *)
h[x_]:=g[2μX-x]; (* define atithetic variate h *)
GG=Expectation[g[x],x\[Distributed]UniformDistribution[{1,8}]]; (* 𝔼[g(X)] *)
HH=Expectation[h[x],x\[Distributed]UniformDistribution[{1,8}]]; (* 𝔼[h(x)] *)
varg=Expectation[((g[x])^2),x\[Distributed]UniformDistribution[{1,8}]]-(GG^2); (* var(g(X)) *)
varh=Expectation[((h[x])^2),x\[Distributed]UniformDistribution[{1,8}]]-(HH^2); (* var(g(X)) *)
covgh=Expectation[g[x]h[x],x\[Distributed]UniformDistribution[{1,8}]]-GG HH;  (* cov(g(X),h(X)) *)
varGnTilde=((varg+varh+2covgh)/(4n)); (* var((Subscript[(OverTilde[G]), n])) *)
varGn=(varg/(2n)); (* var((Subscript[G, 2n])) *)
(varGnTilde/varGn)//N (* ratio var((Subscript[(OverTilde[G]), n])) to (Subscript[G, n]) *)

Using Mathematica

n=10^5; (* sample size *)
g[x_]:=21x (E^(-x/2)); (* define function g *)
μX=Mean[UniformDistribution[{1,8}]]; (* mean of U(1,8) *)

rv=RandomReal[UniformDistribution[{1,8}],n]; (* generate U(1,8) RVs *)
grv=g[rv]; (* apply g to RVs *)
hrv=g[2 μX-rv]; (* apply g to antithetic RVs *)
rv=(grv+hrv)/2; (* take average *)

GnHat=Mean[rv]; (* calculate mean of estimator *)
varGnHat=Variance[rv]/n;  (* calculate sample variance of estimator *)
sdGnHat=(Sqrt[varGnHat]); (* calculate sample standard error *)
{GnHat,varGnHat,sdGnHat}//N (* display results *)

Using R

See R code file.

(c)

The RV    X~U(1,8)    has characteristic function
(Subscript[φ, X])(u) = 𝔼[(e^(i u X))] = (∫_(-∞)^∞)(e^(i u x))(1/7)I(1≤x≤8)dx
  = (1/7)(∫_1^8)(e^(i u x))dx = (Subscript[(1/7)[(1/(i u))(e^(i u x))], 1]^8)
  = (1/(7i u))((e^(8i u))-(e^(i u)))
when u≠0 (equals 1 when u=0).

The RV    2𝔼[X]-X    has CF
𝔼[(e^(i u(2𝔼[X]-X)))] = (e^(2i u 𝔼[X]))𝔼[(e^(-i u X))]
  = (e^(9i u))(Subscript[φ, X])(-u)
  = -(1/(7i u))(e^(9i u))((e^(-8i u))-(e^(-i u)))
  = (1/(7i u))((e^(8i u))-(e^(i u)))
when u≠0 (equals 1 when u=0).

As the CFs are the same    X    and    2𝔼[X]-X    have the same distribution.

Calculations

Mean[UniformDistribution[{1,8}]]

CharacteristicFunction[UniformDistribution[{1,8}],u]//FullSimplify

Expectation[(E^(I u(9-x))),x\[Distributed]UniformDistribution[{1,8}]]//FullSimplify

Q2

(a)

The distribution function of X is
(Subscript[F, X])(x) = (∫_(-∞)^x)(Subscript[f, X])(x)dx
 = 24(∫_0^x)(((3z+2))^(-3))dz = (Subscript[24[-(1/6)(((3z+2))^(-2))], 0]^x)
  =-(Subscript[4[(((3z+2))^(-2))], 0]^x) = -4((((3x+2))^(-2))-(1/4))
  = 1-(4/(((3x+2))^2)) .

Now set
u = (Subscript[F, X])(x)
 = 1-(4/(((3x+2))^2)) = (((((3x+2))^2)-4)/(((3x+2))^2))
  = ((9(x^2)+12x)/(9(x^2)+12x+4))
or
u(9(x^2)+12x+4) = 9(x^2)+12x
which is
9(u-1)(x^2)+12(u-1)x+4u = 0.

Using the quadratic formula to solve for x we have
x = ((-12(u-1)±(Sqrt[(((12(u-1)))^2)-4*9(u-1)4u]))/(2*9(u-1)))
  = ((-12(u-1)±(Sqrt[144((u^2)-2u+1)-144((u^2)-u)]))/(18(u-1)))
  = ((-12(u-1)±12(Sqrt[1-u]))/(18(u-1))) = ((-2(u-1)±2(Sqrt[1-u]))/(3(u-1)))
  = -(2/3)±(2/3)((Sqrt[1-u])/(u-1))
  = -(2/3)∓(2/(3(Sqrt[1-u])))
for 0<u<1.

However, the condition x≥0 means the only valid solution is
x = -(2/3)+(2/(3(Sqrt[1-u])))
or
(Subscript[F, X]^(-1))(u) = (Subscript[F, X]^(-1))((Subscript[F, X])(x)) = x = (2/(3(Sqrt[1-u])))-(2/3).

Calculations

fX[x_]:=(24/(((3x+2))^3));
Assuming[{x>0},(∫_0^∞)fX[x]\[DifferentialD]x] (* check (Subscript[f, X]) integrates to 1 *)
Assuming[{x>0},(∫_0^x)fX[z]\[DifferentialD]z]//FullSimplify  (* CDF (Subscript[F, X]) *)
(∫_0^∞)x*fX[x]\[DifferentialD]x//FullSimplify  (* 𝔼[X] *)

Roots[9(u-1)(x^2)+12(u-1)x+4u==0,x]//FullSimplify

Plot[-(2/3)-(2/(3(Sqrt[1-u]))),{u,0,1},AxesLabel->{"u","x"}]

Plot[-(2/3)+(2/(3(Sqrt[1-u]))),{u,0,1},AxesLabel->{"u","x"}]

(b)

If U~Uniform(0,1) we have by method page 41 Chapter 3 Notes
(Subscript[F, X]^(-1))(U) = (2/(3(Sqrt[1-U])))-(2/3) ~ X.

Using Mathematica

n=(10^5) ; (* sample size *)
fX[x_]:=(24/(((3x+2))^3)); (* PDF *)
fInv[u_]:=(2/(3(Sqrt[1-u])))-(2/3); (* inverse CDF *)
rv=RandomReal[UniformDistribution[],n]; (* generate U(0,1) RVs *)
rv=Map[fInv,rv]; (* apply inverse CDF *)
histogram=Histogram[rv,Automatic,"PDF"];
density=Plot[fX[x],{x,0,4},PlotRange->All,PlotStyle->{Blue,Thick}];
Print[Show[{histogram,density}]];
{Max[rv],Mean[rv],StandardDeviation[rv]/(Sqrt[n])}

Using R

See R code file.

Q3

(a)

We have is Gaussian with
𝔼[(Subscript[X, t])] = -(t^2)
and
Q(s,t) = cov((Subscript[X, s]),(Subscript[X, t]))=(e^(-4|t-s|))cos((π/4)|t-s|)).

(Subscript[X, t]) is Guassian with
𝔼[(Subscript[X, 2])] = -4
and
var((Subscript[X, 2])) = Q(2,2) = 1
so the density
(Subscript[f, (Subscript[X, 3])])(x) = (1/(Sqrt[2π Q(2,2)]))exp(-((((x-𝔼[(Subscript[X, 2])]))^2)/(2Q(2,2))))
  = (1/(Sqrt[2π]))exp(-((((x+4))^2)/2))
and characteristic function
(Subscript[φ, (Subscript[X, 7])])(u) = exp(i u 𝔼[(Subscript[X, 2])] - (1/2)(Q^2)(2,2)(u^2))
  = exp(-4i u - (1/2)(u^2)).

Calculations

PDF[NormalDistribution[-4,1],x]
CharacteristicFunction[NormalDistribution[-4,1],u]

(b)

By TNC Part 3 (page 19 Chapter 2 Notes)
𝔼[(({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))|(Subscript[X, 1])] = 𝔼[(({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))]+cov((({{(Subscript[X, 3])}, {(Subscript[X, 5])}})),(Subscript[X, 1]))cov((((Subscript[X, 1]),(Subscript[X, 1])))^(-1))(({{(Subscript[X, 1])-𝔼[(Subscript[X, 1])]}}))
  = (({{𝔼[(Subscript[X, 3])]}, {𝔼[(Subscript[X, 5])]}}))+(({{cov((Subscript[X, 3]),(Subscript[X, 1]))}, {cov((Subscript[X, 5]),(Subscript[X, 1]))}}))var((((Subscript[X, 1])))^(-1))((Subscript[X, 1])-𝔼[(Subscript[X, 1])])
  = (({{-9}, {-25}}))+(({{-(e^(-8))}, {(e^(-16))}}))*1*((Subscript[X, 1])+1)
  = (({{-(e^(-8))}, {(e^(-16))}}))(Subscript[X, 1])+(({{-9-(e^(-8))}, {-25+(e^(-16))}})).

We can also use scalar form of TNC as
𝔼[(({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))|(Subscript[X, 1])] = (({{𝔼[(Subscript[X, 3])|(Subscript[X, 1])]}, {𝔼[(Subscript[X, 3])|(Subscript[X, 1])]}})) 
where
𝔼[(Subscript[X, 3])|(Subscript[X, 1])] = 𝔼[(Subscript[X, 3])]+((cov((Subscript[X, 3]),(Subscript[X, 1])))/(var((Subscript[X, 1]))))((Subscript[X, 1])-𝔼[(Subscript[X, 1])])
  = -9+((-(e^(-8)))/1)((Subscript[X, 1])+1)
  = -(e^(-8))(Subscript[X, 1])-9-(e^(-8))
and
𝔼[(Subscript[X, 5])|(Subscript[X, 1])] = 𝔼[(Subscript[X, 5])]+((cov((Subscript[X, 5]),(Subscript[X, 1])))/(var((Subscript[X, `]))))((Subscript[X, 1])-𝔼[(Subscript[X, 1])])
  = -25+((e^(-16))/1)((Subscript[X, 1])+1)
  = (e^(-16))(Subscript[X, 1])-25+(e^(-16)).

Calculations

cov[s_,t_]:=(E^(-4Abs[t-s]))Cos[(π/2)Abs[t-s]];
{cov[3,1],cov[5,1],cov[1,1]}

(c)

By TNC Part 4 (page 19 Chapter 2 Notes)
var((({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))|(Subscript[X, 1])]) = 𝔼[((({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))-𝔼[(({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))|(Subscript[X, 1])])((((({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))-𝔼[(({{(Subscript[X, 3])}, {(Subscript[X, 5])}}))|(Subscript[X, 1])]))^T)|(Subscript[X, 1])]
  = cov((({{(Subscript[X, 3])}, {(Subscript[X, 5])}})),(({{(Subscript[X, 3])}, {(Subscript[X, 5])}})))-cov((({{(Subscript[X, 3])}, {(Subscript[X, 5])}})),(Subscript[X, 1]))cov((((Subscript[X, 1]),(Subscript[X, 1])))^(-1))cov((((({{(Subscript[X, 3])}, {(Subscript[X, 5])}})),(Subscript[X, 1])))^T)
  = (({{cov((Subscript[X, 3]),(Subscript[X, 3])), cov((Subscript[X, 3]),(Subscript[X, 5]))}, {cov((Subscript[X, 5]),(Subscript[X, 3])), cov((Subscript[X, 5]),(Subscript[X, 5]))}}))-(({{cov((Subscript[X, 3]),(Subscript[X, 1]))}, {cov((Subscript[X, 5]),(Subscript[X, 1]))}}))var((((Subscript[X, 1])))^(-1))(((({{cov((Subscript[X, 3]),(Subscript[X, 1]))}, {cov((Subscript[X, 5]),(Subscript[X, 1]))}})))^T)
  = (({{var((Subscript[X, 3])), cov((Subscript[X, 3]),(Subscript[X, 5]))}, {cov((Subscript[X, 5]),(Subscript[X, 3])), var((Subscript[X, 5]))}}))-(({{cov((Subscript[X, 3]),(Subscript[X, 1]))}, {cov((Subscript[X, 5]),(Subscript[X, 1]))}}))var((((Subscript[X, 1])))^(-1))(({{cov((Subscript[X, 3]),(Subscript[X, 1])), cov((Subscript[X, 5]),(Subscript[X, 1]))}}))
  = (({{1, -(e^(-8))}, {-(e^(-8)), 1}}))-(({{-(e^(-8))}, {(e^(-16))}}))*1*(({{-(e^(-8)), (e^(-16))}}))
  = (({{1, -(e^(-8))}, {-(e^(-8)), 1}}))- (({{(e^(-16)), -(e^(-24))}, {-(e^(-24)), (e^(-32))}}))
  = (({{1-(e^(-16)), -(e^(-8))+(e^(-24))}, {-(e^(-8))+(e^(-24)), 1-(e^(-32))}})).

Calculations

{cov[3,5],cov[3,1],cov[5,1]}

(({{1, -(E^(-8))}, {-(E^(-8)), 1}}))-(({{-(E^(-8))}, {(E^(-16))}})).(({{-(E^(-8)), (E^(-16))}}))//FullSimplify//MatrixForm