RATS 10.1
RATS 10.1

# EMEXAMPLE.RPF

EMEXAMPLE.RPF demonstrates the EM optimization algorithm, applied to a "mixture" model. From the User's Guide, Section 4.14.

Full Program

all 6237
open data g10xrate.xls
data(format=xls,org=columns) / usxjpn usxfra usxsui
*
* This estimates a "mixture" model, with a zero mean and two branches
* for the variance.
*
set xjpn = 100.0*log(usxjpn/usxjpn{1})
compute rstart=2,rend=6237
*
* Get guess values based on the sample variance.
*
stats xjpn
compute lowv=.5*%variance,highv=2.0*%variance
******************************************************************
*
* Direct maximization
*
nonlin p sigsq1 sigsq2
frml logl = log(p*exp(%logdensity(sigsq1,xjpn))+$(1-p)*exp(%logdensity(sigsq2,xjpn))) * compute t1=%cputime() compute p=.5 compute sigsq1=lowv,sigsq2=highv maximize logl rstart rend ****************************************************************** * * Classic "EM" calculation * compute t2=%cputime() compute p=.5 compute sigsq1=lowv,sigsq2=highv compute oldparms=||p,sigsq1,sigsq2|| do iters=1,500 set pstar = f1=exp(%logdensity(sigsq1,xjpn)),$
f2=exp(%logdensity(sigsq2,xjpn)),$p*f1/(p*f1+(1-p)*f2) sstats(mean) rstart rend pstar>>p pstar*xjpn^2>>sumsq1$
(1-pstar)*xjpn^2>>sumsq2
compute sigsq1=sumsq1/p,sigsq2=sumsq2/(1-p)
compute newparms=||p,sigsq1,sigsq2||
if %testdiff(newparms,oldparms)<1.e-6
break
compute oldparms=newparms
end
disp "Estimates from classical EM" p sigsq1 sigsq2
disp "Iterations" iters
******************************************************************
*
* EM using the ESTEP option
*
function EStep
set pstar = f1=exp(%logdensity(sigsq1,xjpn)),$f2=exp(%logdensity(sigsq2,xjpn)),$
p*f1/(p*f1+(1-p)*f2)
end EStep
*
compute t3=%cputime()
frml emlogl = pstar*%logdensity(sigsq1,xjpn)+$(1-pstar)*%logdensity(sigsq2,xjpn)+$
pstar*log(p)+(1-pstar)*log(1-p)
*
compute p=.5
compute sigsq1=lowv,sigsq2=highv
nonlin p sigsq1 sigsq2
maximize(estep=estep(),method=bhhh,cvcrit=1.e-6,iters=500) emlogl rstart rend
compute t4=%cputime()
*
disp "Timing"
disp "ML" @15 t2-t1
disp "Classical EM" @15 t3-t2
disp "Max(ESTEP)" @15 t4-t3

Output

This does maximum likelihood, plus EM two ways (one explicitly handling the E step, the other using the ESTEP option on MAXIMIZE). All three methods come up with the same results, but maximum likelihood, which is easiest to set up, is also the fastest to execute for a model this simple.

Statistics on Series XJPN

Observations                  6236

Sample Mean               0.014351      Variance                   0.398893

Standard Error            0.631580      SE of Sample Mean          0.007998

t-Statistic (Mean=0)      1.794334      Signif Level (Mean=0)      0.072808

Skewness                  0.784016      Signif Level (Sk=0)        0.000000

Kurtosis (excess)        13.265486      Signif Level (Ku=0)        0.000000

Jarque-Bera           46362.538224      Signif Level (JB=0)        0.000000

MAXIMIZE - Estimation by BFGS

Convergence in    15 Iterations. Final criterion was  0.0000004 <=  0.0000100

Usable Observations                      6236

Function Value                     -5316.1941

Variable                        Coeff      Std Error      T-Stat      Signif

************************************************************************************

1.  P                            0.6435210012 0.0252360550     25.50006  0.00000000

2.  SIGSQ1                       0.0990984333 0.0075984972     13.04185  0.00000000

3.  SIGSQ2                       0.9404851069 0.0517637914     18.16878  0.00000000

Estimates from classical EM       0.64351       0.09909       0.94046

Iterations 217

MAXIMIZE - Estimation by BHHH

Convergence in   169 Iterations. Final criterion was  0.0000009 <=  0.0000010

Usable Observations                      6236

Function Value                     -8203.9312

Variable                        Coeff      Std Error      T-Stat      Signif

************************************************************************************

1.  P                            0.6435127797 0.0135557635     47.47153  0.00000000

2.  SIGSQ1                       0.0990960844 0.0040509789     24.46226  0.00000000

3.  SIGSQ2                       0.9404679960 0.0162516182     57.86919  0.00000000

Timing

ML            0

Classical EM  2

Max(ESTEP)    6