MA(1), n=50

A.I. McLeod
Department of Statistical and Actuarial Sciences,
University of Western Ontario,
e-mail: aim@uwo.ca

B. Quenneville,
Statistics Canada
e-mail: quenne@statcan.ca

MA(1) Simulation

These functions are used in the simulation study for the case n=50 as well as in testing the general algorithms for MA(1) mle and mele.

Needs["Statistics`ContinuousDistributions`"];

\(SimulateMA1[n_, \ θ_] := \ 
      With[{a = RandomArray[NormalDistribution[0, 1], \ n + 1]}, 
        Drop[a, {1}] - θ\ Drop[a, {\(-1\)}]]; 

    SimulateMA1::usage\  = 
	
      "SimulateMA1[n,θ] simulates a series of length n from a MA(1) Gaussian process with parameter θ, mean zero and unit innovation variance."; \)

Efficient Numeric Likelihood Computations in Mathematica

Derivation of Algorithms

Newbold's algorithm can be made much more efficient when only numerical values of the log likelihood are needed by using the Mathematica Compiler and by re-writing the calculations involved to make more use of efficient Mathematica functions such as NestList, FoldList and Apply. First consider the computation of the vector L z which is of length n + 1.  After some simplications, we see that \(TraditionalForm\`L\ z\  = \ \(\((\ \ α\_j)\)\^'\)\), where (Subscript[α, 0]) = 0 is the first element and the remaining elements are defined recursively by \(TraditionalForm
      \`α\_j\  = θ\ α\_\(j - 1\)\  + \ Z\_j, \ j = 1, \ 
      2, \  ... \ , \ n\), where (Subscript[Z, 0])  = 0.  This computation is efficiently performed by Mathematica 's FoldList. Since we are at this point we are just interested in numerical evaluation we use the compile function.

GetLz = 
      Compile[{{θ, _Real}, {z,  _Real,  1}}, 
        FoldList[(#1 θ  +  #2)&, 0, z]];

The determinant, \(TraditionalForm
      \`D\  = \ \(1 + θ\^2\  + θ\^4 + \  ... \)\  + 
          θ\^\(2  n\)\), is efficiently computed using NestList to generate the individual terms and then summing.

DetMA  = 
      Compile[{{θ, _Real}, {n,  _Integer}}, 
        Apply[Plus, NestList[#1 θ &, 1, n]^2]];

Next, we evaluate the term h L z/D.  Since \(TraditionalForm
      \`h\ L\ z\  = \ 
        \(θ\ α\_1 + \(θ\^2\) α\_2 + \  ... \)\  + 
          \(θ\^n\) α\_n\) we can use Horner's Rule to efficiently compute this sum. Horner's Rule is implemented in Mathematica using the function Fold.

Getu0  = 
      Compile[{{θ, _Real}, {Lz,  _Real,  1}, {detma,  _Real}}, 
        -Fold[#1 θ  +  #2&, 0, Reverse[Lz]]/detma];

The computation of the sum of squares function \(TraditionalForm
      \`\(S \((θ)\)\  = \ 
        \(\((L\ z\  - \ h\ \(h\^'\)\ L\ z/D)\)\^'\)\ 
          \((L\ z\  - \ h\ \(h\^'\)\ L\ z/D)\)\ \)\)is straightforward. We use the Mathematica compiler to optimize the vector computations.

GetSumSq  =  
      Compile[{{θ, _Real}, {Lz,  _Real,  1}, {u,  _Real}, {n, 
             _Integer}}, 
        Apply[Plus, (Lz + NestList[#1 θ &, 1, n] u)^2]];

Finally, the concentrated loglikelihood function is defined. The computation speed is increased by about a factor of 50 times when n = 50 and is even larger for larger n.

\(\(logLMA1F[θ_, \ z_]\  := \ 
      Module[{n = Length[z]}, Lz = GetLz[θ, z]; 
        detma = DetMA[θ, n]; u = Getu0[θ, Lz, detma]; 
		
        S = GetSumSq[θ, Lz, u, n]; 
        \(-\((1/2)\)\)\ Log[detma] - \ \((n/2)\) Log[S/n]]; \)\)

The mean likelihood estimate \(TraditionalForm\`\(θ\&_\ \)\) can be evaluated using NIntegrate.

Meanle[z_] := 
      NIntegrate[θ E^logLMA1F[θ,  z], {θ, -1, 1}]/
        NIntegrate[E^logLMA1F[θ,  z], {θ, -1, 1}]

Notice that in the above expression the loglikelihood function is evaluated separately in both the numerator and denominator. We can save function evaluations by using our own numerical quadrature routine.

SimpsonQuadratureWeights[k_, a_,  b_] := 
      With[{h = (2 k)/3}, {a + (b - a) Range[0, 2 k]/(2  k), 
          Prepend[Append[Drop[Flatten[Table[{4, 2}, {k}]], {-1}], 1], 1]}]

{X, W} = SimpsonQuadratureWeights[100, -1, 1];

\(\(GETMEANLEF = 
      Compile[{{z, \ _Real, \ 1}, {W, \ _Real, \ 1}, {X, \ _Real, \ 1}, {f, 
            \ _Real, \ 1}}, Plus@@\((X\ f)\)/Plus@@f]				\)\)

\(\(\(MEANLEF[z_] := 
      With[{f = Plus@@W\ E^\((\(logLMA1F[#1, z]&\)/@X)\)}, 
        GETMEANLEF[z, W, X, f]]; \)					\)\)

Our tests indicate acceptable accuracy and about a 70% improvement in speed as compared with Mathematica's more sophisticated NIntegrate function.

Bayesian Analysis

As previously noted the mean likelihood is just the usual Bayes estimator with a uniform prior.  However the standard Bayesian analysis (Box and Jenkins, 1976, p.250-258) recommends the use of a Jeffrey's non-informative prior. In the case of the MA(1) model this prior is given by \(TraditionalForm
      \`π(θ) = 1/√\((1 - θ\^2)\)\).  The following function computes the Bayes estimate under a Jeffrey's prior,

Since the prior is infinite at the end-points we have to change the integration limits for our implementation of Simpson's rule

{X2, W2} = SimpsonQuadratureWeights[100, -0.999, 0.999];

\(\(\(BayesMA1[z_] := 
      With[{f = 
            Chop[Plus@@W2\ *\((\(1/Sqrt[1 - #1^2]&\)/@X2)\)*
                E^\((\(logLMA1F[#1, z]&\)/@X2)\)]}, 
        GETMEANLEF[z, W2, X2, f]]; \)					\)\)

SeedRandom[1]; z = SimulateMA1[50,  0.5];

BayesMA1[z]

SeedRandom[1]; Table[BayesMA1[SimulateMA1[50, 0.95]], {30}]

As check on our estimator we can also compute it the long way,

Bayes2[z_] := 
      NIntegrate[
          θ (E^logLMA1F[θ,  z])/Sqrt[1 - (θ^2)], {
            θ, -1, 1}]/
        NIntegrate[
          (E^logLMA1F[θ,  z])/Sqrt[1 - (θ^2)], {θ, 
            -1, 1}]

SeedRandom[1]; Table[Bayes2[SimulateMA1[50, 0.95]], {30}]

Comparison of Methods

First, let's compare results using the two methods of evaluating the general likelihood for the MA(1) model. We use the same random numbers so we can compare the results.

Evaluation of Likelihood

SeedRandom[1]; 
    Timing[Table[LogLikelihoodMA1[0.5, SimulateMA1[50, 0.5]], {10}]]

SeedRandom[1]; 
    Timing[Table[logLMA1F[0.5, SimulateMA1[50, 0.5]], {10}]]

Optimization of Likelihood Function

The elaborate code in MLE is needed to handle the case where the mle is exactly 1 or -1.  In this case, FindMinimum fails and we decide then whether the MLE is 1 or -1.

\(MLE[z_] := 
      Module[{ans}, 
        ans = FindMinimum[
            \(-logLMA1F[theta, z]\), {theta, {\(-0.5\), 0.5}, \(-1\), 1}]; 

        	If[ListQ[ans], 
          theta /. 
            ans[[2]], {\(-1\), 1}
              [[\(With[{
                    t = \(logLMA1F[#1, z]&\)/@{\(-1\), 1}}, 
                  Position[t, Max[t]]] // Flatten\) // First
              ]]]]\)

It is helpful to turn off error messages:

Off[FindMinimum::regex]

SeedRandom[1]; 
    MLE[SimulateMA1[50,  1]] 
      Timing[Table[MLE[SimulateMA1[50,  0.75]], {10}]]

So here is our MA(1) likelihood optimization function,

\(\(MLEMA1[z_]\  := 
      Module[{ETOL = 0.0001, 
          NULLVEC = {{NULL, NULL}, {NULL, NULL}, {NULL, NULL}}}, 
		
        ls = Table[{θ, logLMA1F[θ, z]}, {θ, \(-1\), 
              1,  .1}]; 

        		{maxls}\  = 
          Select[ls, \((#1[\([2]\)] == Max[\(Transpose[ls]\)[\([2]\)]])\)&]; 
        
		i = Floor[maxls[\([1]\)]*10 + 11]; 

        			{{a, fa}, {b, fb}, {c, fc}} = 
          Which[i == 1, \ d = \(- .999\); fd = logLMA1F[d, z]; 
            If[\(ls[\([1]\)]\)[\([2]\)] > fd, 
              NULLVEC, {ls[\([1]\)], {d, fd}, ls[\([2]\)]}], 
				
            i == 21, d = 0.999; fd = logLMA1F[d, z]; 
            If[\(ls[\([21]\)]\)[\([2]\)] > fd, 
              NULLVEC, {ls[\([20]\)], {d, fd}, ls[\([21]\)]}], 
				
            \((i > 1)\) && \((i < 21)\), 				

            		{{a, fa}, {b, fb}, {c, fc}} = ls[\([{i - 1, i, i + 1}]\)]]; 

        		If[a === NULL, Which[i == 1, b = \(-1. \), i == 21, b = 1. ], 

          		While[
            \((test1 = \((c - a)\) > ETOL)\)\  || \ 
              \((\((test2 = \((fb - Max[fa, fc])\))\) > ETOL)\), 
		
            fx = logLMA1F[
                x = b + 
                    0.5 \(\(\((b - a)\)\^2\) \((fb - fc)\) - 
                          \(\((b - c)\)\^2\) \((fb - fa)\)\)\/\(\((b - a)\) 
                            \((fb - fc)\) - \((b - c)\) \((fb - fa)\)\), z]; 
            
		{{a, fa}, {b, fb}, {c, fc}} = 
              If[fx > fb, 
                Which[\((x < b)\), {{a, fa}, {x, fx}, {b, fb}}, 
                  \((x > b)\), {{b, fb}, {x, fx}, {c, fc}}]\ , 
	\ 
                If[\((b - a)\) > 0.5\ \((c - a)\), 
                  d = b + 0.38196*\((a - b)\); fd = logLMA1F[d, z]; 
                  If[fd > 
                      fb, {{a, fa}, {d, fd}, {b, fb}}, {{d, fd}, {b, fb}, {c, 
                        fc}}], d = b + 0.38196*\((c - b)\); 
                  fd = logLMA1F[d, z]; 
                  If[fd > 
                      fb, {{b, fb}, {d, fd}, {c, fc}}, {{a, fa}, {b, fb}, {d, 
                        fd}}]]]]]; b]			\)\)

SeedRandom[1]; 
    MLE[SimulateMA1[50,  1]] 
      Timing[Table[MLEMA1[SimulateMA1[50, 0.75]], {10}]]

We ran the simulations using the FindMinimum approach and our MLEMA1 function.  For parameter values near 1 or -1 we got slightly different results.  It seemed that the MLEMA1 were more accurate and gave results closer to those predicted by Cryer and Ledholter for the probability of getting a parameter on the unit boundary.

Mean Likelihood

The mean likelihood estimate \(TraditionalForm\`\(θ\&_\ \)\) can be evaluated using NIntegrate.

Meanle[z_] := 
      NIntegrate[θ E^logLMA1F[θ,  z], {θ, -1, 1}]/
        NIntegrate[E^logLMA1F[θ,  z], {θ, -1, 1}]

Notice that in the above expression the loglikelihood function is evaluated separately in both the numerator and denominator. We can save function evaluations by using our own numerical quadrature routine.

SimpsonQuadratureWeights[k_, a_,  b_] := 
      With[{h = (2 k)/3}, {a + (b - a) Range[0, 2 k]/(2  k), 
          Prepend[Append[Drop[Flatten[Table[{4, 2}, {k}]], {-1}], 1], 1]}]

{X, W} = SimpsonQuadratureWeights[100, -1, 1];

\(\(GETMEANLEF = 
      Compile[{{z, \ _Real, \ 1}, {W, \ _Real, \ 1}, {X, \ _Real, \ 1}, {f, 
            \ _Real, \ 1}}, Plus@@\((X\ f)\)/Plus@@f]				\)\)

\(\(\(MEANLEF[z_] := 
      With[{f = Plus@@W\ E^\((\(logLMA1F[#1, z]&\)/@X)\)}, 
        GETMEANLEF[z, W, X, f]]; \)					\)\)

Our tests indicate acceptable accuracy and about a 70% improvement in speed as compared with Mathematica's more sophisticated NIntegrate function.

SeedRandom[1]; Timing[Table[Meanle[SimulateMA1[50, 0.5]], {10}]]

SeedRandom[1]; Timing[Table[MEANLEF[SimulateMA1[50, 0.5]], {100}]]

Bayes Estimator

{X2, W2} = SimpsonQuadratureWeights[100, -0.999, 0.999];

\(\(\(BayesMA1[z_] := 
      With[{f = 
            Chop[Plus@@W2\ 
                E^\((\(logLMA1F[#1, z]*1/Sqrt[1 - #1^2]&\)/@X2)\)]}, 
        GETMEANLEF[z, W2, X2, f]]; \)					\)\)

SeedRandom[1]; Timing[Table[BayesMA1[SimulateMA1[50, 0.5]], {100}]]

Variance of the Estimated Efficiency

Suppose N simulations are to be performed and N is quite large. Let (Subscript[X, i]) and (Subscript[Y, i]), i = 1,   ...  ,  N denote the squared error in the i-th simulation for estimators 𝒳 and 𝒴 respectively. We will assume that before the simulations are run, that the bivariate random variables ((Subscript[X, i]),  (Subscript[Y, i ])) are an independent sequence of random variables with finite variance. Note that  (Subscript[X, i]) and (Subscript[Y, i]) may themselves be correlated and hence not independent of each other but that it is the sequence which is assumed to be independent. 
Then  sample average squared errors, \(TraditionalForm\`\((X\&_, \ Y\&_)\)\),  are approximately bivariate normal since by the usual central limit theorem, any linear combination of  \(TraditionalForm\`\(X\&_\ \)\) and  \(TraditionalForm\`Y\&_\) is approximately normal.  The estimated relative efficiency of method 𝒴 vs. method 𝒳 is \(TraditionalForm\`e\&^\  = X\&_/\ Y\&_\).  Let \(μ\_x\&_\)=E\({X\&_}\) and \(μ\_y\&_\)=E\({Y\&_}\). Employing a Taylor series linearization of  \(TraditionalForm\`\(e\&^\ \)\),

\(e\&^\  ≐ 
              μ\_x\&_\/μ\_y\&_ + 
                \((X\&_ - μ\_x\&_)\) 1\/μ\_y\&_ - \ 
                \((\ Y\&_ - μ\_y\&_\ \ )\) μ\_x\&_\/μ\_y\&_\%2\)

Hence, after some algebraic simpliciation it can be shown that,

\(Var {e\&^\ } ≐ 
                \(1\/N\) {
                    σ\_x\/μ\_y\%2 + 
                      \(\(σ\_y\) μ\_x\%2\)\/μ\_y\%\(\ 4\) - \ 
                      2  σ\_\(x, y\)\ \ μ\_x\/μ\_y\%\(\ 3\)\)} ,

where (Subscript[μ, x])=E{(Subscript[X, i])}, (Subscript[μ, y]) =E{(Subscript[Y, i])},  (Subscript[σ, x]) = Var {(Subscript[X, i])}, (Subscript[σ, y]) = Var {(Subscript[Y, i])} and  (Subscript[σ, x,  y]) = Cov {(Subscript[X, i]),  (Subscript[Y, i])}. These values can then be estimated by the sample moments of the observed (Subscript[X, i])'s and (Subscript[Y, i])'s.

The following Mathematica function evaluates the estimated relative efficiency and its standard error for method 𝒴 vs. method 𝒳 given observed data (Subscript[X, i]) and (Subscript[Y, i]), i = 1,   ...  ,  N .

Remove[RelativeEfficiency];

RelativeEfficiency[{X_List, Y_List}] := 
  Module[{μx = Mean[X], μy = Mean[Y], σx, σy, σxy}, 
   σx = Mean[(X - μx)^2]; σy = Mean[(Y - μy)^2]; 
    σxy = Mean[(X - μx)*(Y - μy)]; e = μx/μy; 
    se = Sqrt[(σx/μy^2 + σy*μx^2/μy^4 - 
         2*σxy*μx/μy^3)/Length[X]]; {e, se}];

Simulation Results

Tabulations: MELE vs MLE

Needs["Statistics`Master`"]

RelativeEfficiency[{X_List, Y_List}] := 
  Module[{μx = Mean[X], μy = Mean[Y], σx, σy, σxy}, 
   σx = Mean[(X - μx)^2]; σy = Mean[(Y - μy)^2]; 
    σxy = Mean[(X - μx)*(Y - μy)]; e = μx/μy; 
    se = Sqrt[(σx/μy^2 + σy*μx^2/μy^4 - 
         2*σxy*μx/μy^3)/Length[X]]; {e, se}]; 
  GetSimStudy[base_:""]:=
	Module[{α,i=0,outtitle, outfile,BM=TSUM=R=PMC=Table[{0}, {41}]},
		baseout=StringJoin[base,"out"];
     Do[
       outtitle=StringJoin[baseout,ToString[Chop[10 θ]]];
       outfile=OpenRead[outtitle];
       θs=Transpose[Read[outfile,Expression]]//Transpose;
			Close[outfile];	
			R[[++i]]=RelativeEfficiency[Transpose[(θs-θ)^2]]//N;
						 α=
        N[Length[
              Select[θs, 
                Abs[#1[[1]]-θ]> Abs[#1[[2]]-θ] &]]/
            Length[θs]];
		  PMC[[i]]={α,Sqrt[α (1-α)/Length[θs]]};
		   BM[[i]]={BiasSummary/@Transpose[θs-θ],
          MseSummary/@Transpose[θs-θ]};
	     Print[outtitle,Mean/@Transpose[θs]];
			      ,{θ,-1.0,1.0,0.05}
   ];{Flatten/@BM,R, PMC}];
BiasSummary[x_]:=Mean[Abs[x]]; MseSummary[x_]:=Mean[x^2];

{BM, R, PMC} = 
      GetSimStudy["/home/lexis/aim/p/math/1999/mele/mele/"];

TableForm[
      Transpose[
        Join[{Table[θ, {θ, -1, 1, 0.05}]}, 
          Join[Transpose[R], Transpose[PMC]]]], 
      TableHeadings -> {
          None, {"θ", "R",  "sd(R)", "PMC", 
            "sd(PMC)"}}]

\(TableForm[
      Transpose[
        Join[{Table[θ, {θ, \(-1\), 1, 0.05}]}, Transpose[BM]]], 
      TableHeadings -> {
          None, {"θ", \*""B(\!\(θ\&^\))"", 
            \*""B(\!\(θ\&_\))"", 
            \*""mse(\!\(θ\&^\))"", 
            \*""mse(\!\(θ\&_\))""}}]\)

Tabulations: Bayes vs MLE

Needs["Statistics`Master`"]

RelativeEfficiency[{X_List, Y_List}] := 
  Module[{μx = Mean[X], μy = Mean[Y], σx, σy, σxy}, 
   σx = Mean[(X - μx)^2]; σy = Mean[(Y - μy)^2]; 
    σxy = Mean[(X - μx)*(Y - μy)]; e = μx/μy; 
    se = Sqrt[(σx/μy^2 + σy*μx^2/μy^4 - 
         2*σxy*μx/μy^3)/Length[X]]; {e, se}]; 
  GetSimStudy[base_:""]:=
	Module[{α,i=0,outtitle, outfile,BM=TSUM=R=PMC=Table[{0}, {41}]},
		baseout=StringJoin[base,"out"];
     Do[
       outtitle=StringJoin[baseout,ToString[Chop[10 θ]]];
       outfile=OpenRead[outtitle];
       θs=Transpose[Read[outfile,Expression]]//Transpose;
			Close[outfile];	
			R[[++i]]=RelativeEfficiency[Transpose[(θs-θ)^2]]//N;
						 α=
        N[Length[
              Select[θs, 
                Abs[#1[[1]]-θ]> Abs[#1[[2]]-θ] &]]/
            Length[θs]];
		  PMC[[i]]={α,Sqrt[α (1-α)/Length[θs]]};
		   BM[[i]]={BiasSummary/@Transpose[θs-θ],
          MseSummary/@Transpose[θs-θ]};
	     Print[outtitle,Mean/@Transpose[θs]];
			      ,{θ,-1.0,1.0,0.05}
   ];{Flatten/@BM,R, PMC}];
BiasSummary[x_]:=Mean[Abs[x]]; MseSummary[x_]:=Mean[x^2];

{BM, R, PMC} = 
      GetSimStudy["/home/lexis/aim/p/math/1999/mele/bayes/"];

TableForm[
      Transpose[
        Join[{Table[θ, {θ, -1, 1, 0.05}]}, 
          Join[Transpose[R], Transpose[PMC]]]], 
      TableHeadings -> {
          None, {"θ", "R",  "sd(R)", "PMC", 
            "sd(PMC)"}}]

\(TableForm[
      Transpose[
        Join[{Table[θ, {θ, \(-1\), 1, 0.05}]}, Transpose[BM]]], 
      TableHeadings -> {
          None, {"θ", \*""B(\!\(θ\&^\))"", 
            \*""B(\!\(θ\&_\))"", 
            \*""mse(\!\(θ\&^\))"", 
            \*""mse(\!\(θ\&_\))""}}]\)

\(\(PlotRES[{R_, SD_}, target_] := 
      Module[{crit = 3.29053\  (*99.99  %  C . I . *) , 
          thetas = Table[θ, \ {θ, \(-1\), 1, 0.05}]}, 
        epi = {Thickness[0.01], 
            RGBColor[0.1, 0.1, 0.1], {
              Line/@Transpose[
                  Transpose/@{{thetas, R - crit*\ SD}, {thetas, 
                        R + crit*\ SD}}]}, {RGBColor[0.8, 0.8, 0.8], 
              Thickness[0.001], Line[{{\(-1\), target}, {1, target}}]}}; 

        				ListPlot[Transpose[{thetas, \ R}], 
          PlotStyle -> {PointSize[0.02]}, 
          PlotRange -> {Min[0, 0.94*\((R - crit*\ SD)\)], 
              Max[1, 1.06*\((R + crit*SD)\)]}, AxesOrigin -> {0, 0}, 
          Frame -> True, Epilog -> epi]]; \)\)

PlotRES[PMC // Transpose, 0.5];

PlotRES[R // Transpose, 1];

R >> "r-out.m"

PMC >> "pmc-out.m"

Plots

SetDirectory["c:\\p\\math\\1999\\mele\\"];

rmele = << "rmele.m"; rbayes = << "rbayes.m";

Needs["Graphics`MultipleListPlot`"];

Determine the quantile for a 99.9% confidence interval.

Needs["Statistics`NormalDistribution`"];

Quantile[NormalDistribution[0, 1], 0.995]

\(PlotEFF[rmele_, \ rbayes_] := 
      Module[{x, y, RMELE, SDRMELE, RBAYES, SDBAYES, t, 
          epi = {RGBColor[0.8, 0.8, 0.8], Thickness[0.001], 
              Line[{{\(-1\), 1}, {1, 1}}]}}, 
        t = Table[θ, {θ, \(-1\), 1, 0.05}]; 

        		{RMELE, SDRMELE} = Transpose[rmele]; 

        		{RBAYES, SDBAYES} = Transpose[rbayes]; 
		
        x = MapThread[
              Append, {Transpose[{t, RMELE}], SDRMELE}] /. {a_?NumberQ, 
                b_?NumberQ, c_?NumberQ} -> {{a, b}, ErrorBar[2.57583*c]}; 
	
        y = MapThread[
              Append, {Transpose[{t, RBAYES}], SDBAYES}] /. {a_?NumberQ, 
                b_?NumberQ, c_?NumberQ} -> {{a, b}, ErrorBar[2.57583*c]}; 	

        		MultipleListPlot[x, y, Frame -> True, Axes -> False, 
          PlotRange -> {0.1, 1.9}, 
          SymbolStyle -> {PointSize[0.01], PointSize[0.01]}, 
          PlotJoined -> True, PlotLegend -> {"mele", "Bayes"}, 
          Epilog -> epi]]\)

PlotEFF[rmele, rbayes];

pmcmele = << "pmcmele.m"; pmcbayes = << "pmcbayes.m";

\(PlotPMC[pmcmele_, \ pmcbayes_] := 
      Module[{x, y, PMCMELE, SDPMCMELE, PMCBAYES, SDPMCBAYES, t, 
          epi = {RGBColor[0.8, 0.8, 0.8], Thickness[0.001], 
              Line[{{\(-1\), 0.5}, {1, 0.5}}]}}, 
        t = Table[θ, {θ, \(-1\), 1, 0.05}]; 

        		{PMCMELE, SDPMCMELE} = Transpose[pmcmele]; 

        		{PMCBAYES, SDPMCBAYES} = Transpose[pmcbayes]; 
		
        x = MapThread[
              Append, {Transpose[{t, PMCMELE}], SDPMCMELE}] /. {a_?NumberQ, 
                b_?NumberQ, c_?NumberQ} -> {{a, b}, ErrorBar[2.57583*c]}; 
	
        y = MapThread[
              Append, {Transpose[{t, PMCBAYES}], SDPMCBAYES}] /. {a_?NumberQ, 
                b_?NumberQ, c_?NumberQ} -> {{a, b}, ErrorBar[2.57583*c]}; 	

        		MultipleListPlot[x, y, Frame -> True, Axes -> False, 
          PlotRange -> {0. , 1. }, \ 
          SymbolStyle -> {PointSize[0.01], PointSize[0.01]}, 
          PlotJoined -> True, PlotLegend -> {"mele", "Bayes"}, 
          Epilog -> epi]]\)

PlotPMC[pmcmele, pmcbayes];