NMST 434, Exercise session II: Method of Moments 
February 27, 2020

Example 8: MoM in beta distribution

$Assumptions=α>0&&β>0;
M[m_]=Moment[BetaDistribution[α,β],m]; 
τ[α_,β_]={M[1],M[2]}
MoM = Solve[τ[α,β]=={m1,m2},{α,β}] (* MoM estimators based on X, (X^2)*)
J=Inverse[D[τ[α,β],{{α,β}}]]//Simplify (* Jacobian of the inverse transform *)
Σ=Table[M[i+j]-M[i]M[j],{i,1,2},{j,1,2}]//Simplify (* Variance matrix of (X, (X^2)) *)
J.Σ.Transpose[J]//MatrixForm//Simplify (* asymptotic variance matrix n * AVar((OverHat[α]), (OverHat[β])) *)
$Assumptions=True;

For comparison, without the inverse function theorem.

J2 = D[{α,β}/.MoM,{{m1,m2}}] (* Jacobian of the transform *)
J2/.{m1->M[1],m2->M[2]}//Simplify (* Jacobian of the transform at the first two moments *)

Example: MoM in Gamma distribution

$Assumptions=k>0&&θ>0;
M[m_]=Moment[GammaDistribution[k,θ],m]; 
τ[k_,θ_]={M[1],M[2]}
MoM = Solve[τ[k,θ]=={m1,m2},{k,θ}] (* MoM estimators based on X, (X^2)*)
J=Inverse[D[τ[k,θ],{{k,θ}}]]//Simplify (* Jacobian of the inverse transform *)
Σ=Table[M[i+j]-M[i]M[j],{i,1,2},{j,1,2}]//Simplify (* Variance matrix of (X, (X^2)) *)
J.Σ.Transpose[J]//MatrixForm//Simplify (* asymptotic variance matrix n * AVar((OverHat[k]), (OverHat[θ])) *)

For comparison, computation of the asymptotic variance matrix without the inverse function theorem.

J2 = D[{k,θ}/.MoM[[1]],{{m1,m2}}]/.{m1->M[1],m2->M[2]};
J2.Σ.Transpose[J2]//Simplify//MatrixForm
$Assumptions=True;

Example 15: MoM estimation in Pareto distribution

$Assumptions=α>0&&β>4;
M[m_]=Moment[ParetoDistribution[α,β],m]; 
τ[α_,β_]={M[1],M[2]}
Solve[{((α β)/(-1+β)),(((α^2) β)/(-2+β))}=={m1,m2},{α,β}] //Simplify(* MoM estimators based on X, (X^2)*)
J=Inverse[D[τ[α,β],{{α,β}}]]//Simplify (* Jacobian of the inverse transform *)
Σ=Table[M[i+j]-M[i]M[j],{i,1,2},{j,1,2}]//Simplify (* Variance matrix of (X, (X^2)) *)
J.Σ.Transpose[J]//MatrixForm//Simplify (* asymptotic variance matrix n * AVar((OverHat[α]), (OverHat[β])) *)

Maximum likelihood estimation

data=RandomVariate[ParetoDistribution[1,4],20];
logL[α_,β_]=Total[Log[PDF[ParetoDistribution[α,β],data]]];
Plot3D[logL[α,β],{α,0,Min[data]+1},{β,0,5},AxesLabel->Automatic]
$Assumptions=True;

Example: MoM in skew normal distribution

Skew normal distribution SN(μ,σ,α), where μ∈R is the location parameter, σ∈(0,∞) is the scale parameter, and α∈R is the asymmetry parameter.

$Assumptions=σ>0;
Plot[Table[PDF[SkewNormalDistribution[0,2,α],x],{α,{-3,0.5,4}}]//Evaluate,{x,-5,5},Filling->Axis]
Manipulate[Plot[PDF[SkewNormalDistribution[μ,σ,α],x],{x,-5,5},PlotRange->{0,1},PlotLabel->"Density of the Skew-Normal Distribution"],{{μ,0},-5,5,Appearance->"Labeled"},{{σ,1},0.01,5,Appearance->"Labeled"},{{α,0},-10,10,Appearance->"Labeled"}]
M[k_]=Moment[SkewNormalDistribution[μ,σ,α],k];
CM[k_]=CentralMoment[SkewNormalDistribution[μ,σ,α],k];
Table[M[k],{k,1,3}]
Table[CM[k],{k,1,3}]

Direct method of moments estimators, where E X, E (X^2), and E (X^3) are used.

Solve[{M[1]==m1,M[2]==m2,M[3]==m3},{μ,σ,α}]

Method of moments with E X, Var X, and the skewness Skew X.

((CM[3])/((CM[2])^(3/2)))//Simplify
Solve[{M[1]==m1,CM[2]==m2,((CM[3])/((CM[2])^(3/2)))==m3},{μ,σ,α}]

Semi-Manual computation based on E X, Var X, and Skew X. The skewness depends only on α.

((CM[3])/((CM[2])^(3/2)))//Simplify
Solve[%==m3,α]

Manual computation gives closed form expressions for the estimators of the parameters.

Asymptotic distribution of (OverHat[μ]), ((((OverHat[σ])))^2), (OverHat[α]).

Direct approach. A two-step procedure. 
1. Find the joint asymptotic distribution of the sample mean, sample variance, and the sample skewness. (3-dimensional Δ-theorem).
2. Use the 3-dimensional Δ-theorem again to find the asymptotic distribution of the estimators.

Step 1.

Σ=Table[M[i+j]-M[i]M[j],{i,1,3},{j,1,3}]; (* Variance matrix of (EX, (EX^2), (EX^3)) *)
g[a_,b_,c_]={a,b-(a^2),((c-3 a b+2(a^3))/(((b-(a^2)))^(3/2)))};  (* Function that takes the sample averages of X, (X^2), (X^3) to the sample mean, sample variance, and sample skewness *)
J1=D[g[a,b,c],{{a,b,c}}]/.{a->M[1],b->M[2],c-> M[3]}//Simplify;
V0 = J1.Σ.Transpose[J1]//Simplify; 
MatrixForm[V0](* n AVar(sample mean, sample variance, sample skewness) *)

Check: what does the computation give for α=0? (This is the normal distribution)

V0/.α->0//MatrixForm

Step 2.

δ[a_,b_,c_]=((c^(1/3))/(√((2/π)((c^(2/3))+(((((4-π)/2)))^(2/3)))))); 
h[a_,b_,c_]={a-√((2/π)) δ[a,b,c]√((b/(1-(2/π)((δ[a,b,c])^2)))),(b/(1-(2/π)((δ[a,b,c])^2))),((δ[a,b,c])/(√(1-((δ[a,b,c])^2))))};(* Function that takes the sample mean, variance, and skewness to the estimators *)
J2 =D[h[a,b,c],{{a,b,c}}]/.{a->M[1],b->CM[2],c->((CM[3])/((CM[2])^(3/2)))}//Simplify; (* The Jacobian of function h at the population mean, variance, and skewness *)
V =J2.V0.Transpose[J2]//Simplify;
V//MatrixForm (* Finally, n * the asymptotic variance matrix of our estimators *)
V

Application to the simulated data from the R script.

μR= -5.0063024817859105;
σR=3.0325019891716396;
αR=-1.9521085711315336; (* Point estimators obtained from R *)
V/.{μ-> μR, σ->σR, α->αR}//MatrixForm
V/.{μ-> μR, σ->σR, α->-αR}//MatrixForm (* The result appears to work properly only with α positive, because of (c^(1/3)) in the definition of δ *)
(c^(1/3))/. c->αR

A better approach. Use the inverse function theorem (Theorem 4).
Why do we substitute σ twice in the first line of the following code?

J2b = Inverse[D[{M[1],CM[2],((CM[3])/((CM[2])^(3/2)))}/.σ->√s2,{{μ,s2,α}}]/.s2->(σ^2)]//Simplify;
J2b/.{μ-> μR, σ->σR, α->αR}//N
J2/.{μ-> μR, σ->σR, α->αR}//N
Vfinal = J2b.V0.Transpose[J2b]//Simplify;
Vfinal//MatrixForm
Vfinal/.{μ-> μR, σ->σR, α->αR}//MatrixForm//N

We get (pretty much) the same result in a much more compact form.

Real data analysis - Barolo dataset.

data= Import["C:\\Users\\snagy\\Documents\\R\\Barolo.txt","Data"]; (* Replace the path to where the Barolo.txt file is in your computer *)
data = Flatten[data];
Length[data]
data

Moment[data,1]
CentralMoment[data,2]
((CentralMoment[data,3])/((CentralMoment[data,2])^(3/2)))
MoM = h[Moment[data,1],CentralMoment[data,2],((CentralMoment[data,3])/((CentralMoment[data,2])^(3/2)))](* MoM estimators of the parameters *)
VBarolo = V/.{μ->MoM[[1]],σ -> √MoM[[2]],α->MoM[[3]]} (* n * AVar matrix of the estimators *)
VBarolo2 = Vfinal/.{μ->MoM[[1]],σ -> √MoM[[2]],α->MoM[[3]]} (* n * AVar matrix of the estimators *)
MoM[[3]]+{-1,1}Quantile[NormalDistribution[],.975]√(VBarolo[[3,3]]/Length[data])
(* Asymptotic confidence interval for α *)

Numerical solution to the same problem. 
The asymptotic variance matrix was in the approach using the inverse function theorem obtained without any numerical computations - the only numerical computation needed here is to evaluate the point estimates.

NSolve[{M[1]==Moment[data,1],CM[2]==CentralMoment[data,2],((CM[3])/((CM[2])^(3/2)))==((CentralMoment[data,3])/((CentralMoment[data,2])^(3/2)))},{μ,σ,α},Reals]
(σ^2)/.%[[2]]

Maximum likelihood estimation

Closed form expression for MLE cannot be obtained. Numerical solution.

L[μ_,σ_,α_]=Times @@ PDF[SkewNormalDistribution[μ,σ,α],data];
NMaximize[{L[μ,σ,α],σ>0},{μ,σ,α}]
logL[μ_,σ_,α_]=Total[Log[PDF[SkewNormalDistribution[μ,σ,α],data]]];
MLE = NMaximize[{logL[μ,σ,α],σ>0},{μ,σ,α}]

Observed Fisher information matrix

D2logL = -((D[logL[μ,σ,α],{{μ,σ,α},2}])/(Length[data]))/.MLE[[2]]
Inverse[D2logL]

(α/.MLE[[2]])+{-1,1}Quantile[NormalDistribution[],.975]√(Inverse[D2logL][[3,3]]/Length[data])
(* Asymptotic confidence interval for α *)