MLE in monte carlo

[tommyreeson]

Hi,

Currently I am doing a monte carlo project on hetroskedascity and I am not sure how to go about MLE estimation in the programming.

I have already generated the model just unsure how one would go MLE, specifically generating a logL object and implementing the model into the code.

Any help would be fantastic :)

Thanks

Tommy

Edit: This is the code I have now, forgive the comments I'm working off my lecturers

' OLS Monte Carlo Simulation
!n=30 'number of observations
'set the number of simulations !s in the MC experiment
'or we can call !s the number of artificial data sets
!s=1e2
'any scalars starting with "!" will not be saved in the workfile
'they are just used during the program running

'define the workfile range !r
if !n>!s then
!r=!n
else
!r=!s
endif

%fpath=@linepath
cd %fpath
'cd "H:\mc_exe"

workfile mc u 1 !r
smpl 1 !n
'read the explanatory variable from the file generated by mc_genx.prg into a series called x
'read(skipcol=1,skiprow=3) x10.csv x 'to read only one variable
read(skipcol=1,skiprow=2) x30.csv x x2 'to read two variables

vector(3) b 'create a vector to store the true value of intercept and slope
b(1)=1 'true value of intercept
b(2)=2 'true value of the slope
b(3)=3 'true value of the slope'
scalar sigma=0.1 'true value of sigma
scalar sigmasq=sigma^2 'true value of sigma square

series ey=b(1)+b(2)*x+b(3)*x2 'the expected mean of y


scalar df=!n-3 'degrees of freedom for the t test
scalar a=0.05 'set the level of significance
'calculate the critical values for the t-distribution
scalar tlb=@qtdist(a/2,df) 'the t critical value, 2.5% quantile if a=0.05
scalar tub=@qtdist(1-a/2,df) 'the t critical value, 97.5% quantile if a=0.05

'change the range of the work file and initialize the series to store the estimated values
smpl 1 !s
series b1hat=0 'this is to store the estimated values of intercept
series seb1hat=0 'this is to store the estimated standard error of the intercept estimator
series b2hat=0 'this is to store the estimated values of slope
series seb2hat=0 'this is to store the estimated standard error of the slope estimator
series b3hat=0 'this is to store the estimated values of intercept
series seb3hat=0 'this is to store the estimated standard error of the intercept estimator
series sigmahat=0 'this is to store the estimated values of sigma
series sigmasqhat=0 'this is to store the estimated values of sigma square
scalar t1errb2=0 'this is to store how many Type I errors committed
scalar t1errb3=0 'this is to store how many Type I errors committed
scalar bppower=0
scalar whpower=0

series b1hatwls=0 'this is to store the estimated values of intercept
series seb1hatwls=0 'this is to store the estimated standard error of the intercept estimator
series b2hatwls=0 'this is to store the estimated values of slope
series seb2hatwls=0 'this is to store the estimated standard error of the slope estimator
series b3hatwls=0 'this is to store the estimated values of intercept
series seb3hatwls=0 'this is to store the estimated standard error of the intercept estimator
series sigmahatwls=0 'this is to store the estimated values of sigma
series sigmasqhatwls=0 'this is to store the estimated values of sigma square
scalar t1errb2wls=0 'this is to store how many Type I errors committed
scalar t1errb3wls=0 'this is to store how many Type I errors committed
scalar bppowerwls=0
scalar whpowerwls=0

series b1hath=0 'this is to store the estimated values of intercept
series seb1hath=0 'this is to store the estimated standard error of the intercept estimator
series b2hath=0 'this is to store the estimated values of slope
series seb2hath=0 'this is to store the estimated standard error of the slope estimator
series b3hath=0 'this is to store the estimated values of intercept
series seb3hath=0 'this is to store the estimated standard error of the intercept estimator
series sigmahath=0 'this is to store the estimated values of sigma
series sigmasqhath=0 'this is to store the estimated values of sigma square
scalar t1errb2h=0 'this is to store how many Type I errors committed
scalar t1errb3h=0 'this is to store how many Type I errors committed
scalar bppowerh=0
scalar whpowerh=0

series y=0 'initialize the dependent variable
'the for loop for simulations, do it !s times
for !i=1 to !s
smpl 1 !n
'This section generates stuff for the bptest and white test
y=ey+@nrnd*sigma*(exp(x2))^0.5 'generate the dependent variable
equation eq.ls y c x x2 'create an equation call eq and run the ols ignoring hetroskedascity
series res01=resid
'run the regression needed for the BP test (or LM test)
equation eqbp.ls res01^2 c x x2
scalar bpstat=eqbp.@regobs*eqbp.@r2 'the BP (LM) test statistic
scalar bppval=1-@cchisq(bpstat,2) 'calculate the p-value
if(bppval<a)then
bppower=bppower+1
endif

equation eqwh.ls res01^2 c x x2 x^2 x2^2 x*x2
scalar whstat=eqwh.@regobs*eqwh.@r2 'the white (lm) test stat
scalar whval=1-@cchisq(whstat,2) 'calculate the p-value
if(whval<a)then
whpower=whpower+1
endif


y=ey+@nrnd*sigma*(exp(x2))^0.5 'generate the dependent variable
equation eq03.ls y/(exp(x2)*0.5) 1/(exp(x2)*0.5) x/(exp(x2)*0.5) x2/(exp(x2)*0.5) 'create an equation call eq03 and run the wls
series res01wls=resid
'run the regression needed for the BP test (or LM test)
equation eqbpwls.ls res01wls^2 c x x2
scalar bpstatwls=eqbpwls.@regobs*eqbpwls.@r2 'the BP (LM) test statistic
scalar bppvalwls=1-@cchisq(bpstatwls,2) 'calculate the p-value
if(bppvalwls<a)then
bppowerwls=bppowerwls+1
endif

equation eqwhwls.ls res01wls^2 c x x2 x^2 x2^2 x*x2
scalar whstatwls=eqwhwls.@regobs*eqwhwls.@r2 'the white (lm) test stat
scalar whvalwls=1-@cchisq(whstatwls,2) 'calculate the p-value
if(whvalwls<a)then
whpowerwls=whpowerwls+1
endif

y=ey+@nrnd*sigma*(exp(x2))^0.5 'generate the dependent variable
equation eq02.ls(h) y c x x2 'create an equation call eq and run the ols ignoring hetroskedascity
series res01h=resid
'run the regression needed for the BP test (or LM test)
equation eqbph.ls res01h^2 c x x2
scalar bpstath=eqbph.@regobs*eqbph.@r2 'the BP (LM) test statistic
scalar bppvalh=1-@cchisq(bpstath,2) 'calculate the p-value
if(bppvalh<a)then
bppowerh=bppowerh+1
endif

equation eqwhh.ls res01h^2 c x x2 x^2 x2^2 x*x2
scalar whstath=eqwhh.@regobs*eqwhh.@r2 'the white (lm) test stat
scalar whvalh=1-@cchisq(whstath,2) 'calculate the p-value
if(whvalh<a)then
whpowerh=whpowerh+1
endif

'This section generates values for all of the different models

smpl 1 !n
y=ey+@nrnd*sigma*(exp(x2))^0.5 'generate the dependent variable
equation eq.ls y c x x2 'create an equation call eq and run the ols ignoring hetroskedascity
smpl 1 !s
b1hat(!i) = eq.@coefs(1) 'store the intercept estimate
seb1hat(!i)=eq.@stderrs(1) 'store the intercept estimated standard deviation
b2hat(!i) = eq.@coefs(2) 'store the slope estimate
seb2hat(!i)=eq.@stderrs(2) 'store the slope estimated standard deviation
b3hat(!i) = eq.@coefs(3) 'store the intercept estimate
seb3hat(!i)=eq.@stderrs(3) 'store the intercept estimated standard deviation
sigmahat(!i)=eq.@se 'store the estimated standard deviation of regression
sigmasqhat(!i)=sigmahat(!i)^2 'store the estimated variance of regression



smpl 1 !n
y=ey+@nrnd*sigma*(exp(x2))^0.5 'generate the dependent variable
equation eq03.ls y/(exp(x2)*0.5) 1/(exp(x2)*0.5) x/(exp(x2)*0.5) x2/(exp(x2)*0.5) 'create an equation call eq03 and run the wls
smpl 1 !s
b1hatwls(!i) = eq03.@coefs(1) 'store the intercept estimate
seb1hatwls(!i)=eq03.@stderrs(1) 'store the intercept estimated standard deviation
b2hatwls(!i) = eq03.@coefs(2) 'store the slope estimate
seb2hatwls(!i)=eq03.@stderrs(2) 'store the slope estimated standard deviation
b3hatwls(!i) = eq03.@coefs(3) 'store the intercept estimate
seb3hatwls(!i)=eq03.@stderrs(3) 'store the intercept estimated standard deviation
sigmahatwls(!i)=eq03.@se 'store the estimated standard deviation of regression
sigmasqhatwls(!i)=sigmahatwls(!i)^2 'store the estimated variance of regression

smpl 1 !n
y=ey+@nrnd*sigma*(exp(x2))^0.5 'generate the dependent variable
equation eq02.ls(h) y c x x2 'create an equation call eq and run the ols accounting for hetroskedascity
smpl 1 !s
b1hath(!i) = eq02.@coefs(1) 'store the intercept estimate
seb1hath(!i)=eq02.@stderrs(1) 'store the intercept estimated standard deviation
b2hath(!i) = eq02.@coefs(2) 'store the slope estimate
seb2hath(!i)=eq02.@stderrs(2) 'store the slope estimated standard deviation
b3hath(!i) = eq02.@coefs(3) 'store the intercept estimate
seb3hath(!i)=eq02.@stderrs(3) 'store the intercept estimated standard deviation
sigmahath(!i)=eq02.@se 'store the estimated standard deviation of regression
sigmasqhath(!i)=sigmahath(!i)^2 'store the estimated variance of regression



'This section generates stuff for hypothesis testing

'calculate the lower and upper bound under H0:slope=slope_true given the level of significance
'if the estimated slope is outside the bounds, we reject H0 and make a Type I error
!slb2=b(2)+seb2hat(!i)*tlb 'the lower bound
!sub2=b(2)+seb2hat(!i)*tub 'the upper bound
if (b2hat(!i)<!slb2) or (b2hat(!i)>!sub2) then
t1errb2=t1errb2+1 'if we reject H0:slope=slope_true, we make a Type I error
endif 'if condition finishes here
!slb3=b(3)+seb3hat(!i)*tlb 'the lower bound
!sub3=b(3)+seb3hat(!i)*tub 'the upper bound
if (b3hat(!i)<!slb3) or (b3hat(!i)>!sub3) then
t1errb3=t1errb3+1 'if we reject H0:slope=slope_true, we make a Type I error
endif 'if condition finishes here

'calculate the lower and upper bound under H0:slope=slope_true given the level of significance
'if the estimated slope is outside the bounds, we reject H0 and make a Type I error
!slb2wls=b(2)+seb2hatwls(!i)*tlb 'the lower bound
!sub2wls=b(2)+seb2hatwls(!i)*tub 'the upper bound
if (b2hatwls(!i)<!slb2wls) or (b2hatwls(!i)>!sub2wls) then
t1errb2wls=t1errb2wls+1 'if we reject H0:slope=slope_true, we make a Type I error
endif 'if condition finishes here
!slb3wls=b(3)+seb3hatwls(!i)*tlb 'the lower bound
!sub3wls=b(3)+seb3hatwls(!i)*tub 'the upper bound
if (b3hatwls(!i)<!slb3wls) or (b3hatwls(!i)>!sub3wls) then
t1errb3wls=t1errb3wls+1 'if we reject H0:slope=slope_true, we make a Type I error
endif 'if condition finishes here

'calculate the lower and upper bound under H0:slope=slope_true given the level of significance
'if the estimated slope is outside the bounds, we reject H0 and make a Type I error
!slb2h=b(2)+seb2hath(!i)*tlb 'the lower bound
!sub2h=b(2)+seb2hath(!i)*tub 'the upper bound
if (b2hath(!i)<!slb2h) or (b2hath(!i)>!sub2h) then
t1errb2h=t1errb2h+1 'if we reject H0:slope=slope_true, we make a Type I error
endif 'if condition finishes here
!slb3h=b(3)+seb3hath(!i)*tlb 'the lower bound
!sub3h=b(3)+seb3hath(!i)*tub 'the upper bound
if (b3hath(!i)<!slb3h) or (b3hath(!i)>!sub3h) then
t1errb3h=t1errb3h+1 'if we reject H0:slope=slope_true, we make a Type I error
endif 'if condition finishes here




next 'for loop finishes here
'This section reflects all of our stats from the simulations and puts it into probablity
t1errb2=t1errb2/!s 'percentage of Type I errors committed
t1errb3=t1errb3/!s 'percentage of Type I errors committed
bppower=bppower/!s
whpower=whpower/!s
t1errb2wls=t1errb2wls/!s 'percentage of Type I errors committed in wls
t1errb3wls=t1errb3wls/!s 'percentage of Type I errors committed in wls
bppowerwls=bppowerwls/!s
whpowerwls=whpowerwls/!s
t1errb2h=t1errb2h/!s 'percentage of Type I errors committed when we assume hetroskedascity
t1errb3h=t1errb3h/!s 'percentage of Type I errors committed
bppowerh=bppowerh/!s
whpowerh=whpowerh/!s
'the following is to calculate the power of the t test based on different (wrong) test values
!mintv=0 'minimum test value
!maxtv=4 'maximum test value
!ntestv=101 'number of test values
!step=(!maxtv-!mintv)/!ntestv 'step size
'create a !ntestv by 2 matrix called ptest
'the first column stores the wrong test values while the second stores the power of tests
matrix(!ntestv,2) pcurveb2
matrix(!ntestv,2) pcurveb2wls
matrix(!ntestv,2) pcurveb2h
series reject=0 'initialize a series to store the occurrences of rejecting the wrong hypothesis
series rejectwls=0 'initialize a series to store the occurrences of rejecting the wrong hypothesis wls
series rejecth=0 'initialize a series to store the occurrences of rejecting the wrong hypothesis assuming hetroskedascity

for !i=1 to !ntestv
pcurveb2(!i,1)=!mintv+(!i-1)*!step 'the test value
reject=@recode(((b2hat-pcurveb2(!i,1))/seb2hat<tlb) or ((b2hat-pcurveb2(!i,1))/seb2hat>tub), 1, 0 )
pcurveb2(!i,2)=@sum(reject)/!s
pcurveb2wls(!i,1)=!mintv+(!i-1)*!step 'the test value
rejectwls=@recode(((b2hatwls-pcurveb2wls(!i,1))/seb2hatwls<tlb) or ((b2hatwls-pcurveb2wls(!i,1))/seb2hatwls>tub), 1, 0 )
pcurveb2wls(!i,2)=@sum(rejectwls)/!s
pcurveb2h(!i,1)=!mintv+(!i-1)*!step 'the test value
rejecth=@recode(((b2hath-pcurveb2h(!i,1))/seb2hath<tlb) or ((b2hath-pcurveb2h(!i,1))/seb2hath>tub), 1, 0 )
pcurveb2h(!i,2)=@sum(rejecth)/!s
next
pcurveb2.xy 'plot the power of the t test curve column 2 (c2) against the test values in column 1 (c1)
!mintv=0 'minimum test value
!maxtv=6 'maximum test value
!ntestv=101 'number of test values
!step=(!maxtv-!mintv)/!ntestv 'step size
'create a !ntestv by 2 matrix called ptest
'the first column stores the wrong test values while the second stores the power of tests
matrix(!ntestv,2) pcurveb3
series reject=0 'initialize a series to store the occurrences of rejecting the wrong hypothesis

pcurveb2wls.xy 'plot the power of the t test curve column 2 (c2) against the test values in column 1 (c1)
!mintv=0 'minimum test value
!maxtv=6 'maximum test value
!ntestv=101 'number of test values
!step=(!maxtv-!mintv)/!ntestv 'step size
'create a !ntestv by 2 matrix called ptest
'the first column stores the wrong test values while the second stores the power of tests
matrix(!ntestv,2) pcurveb3wls
series rejectwls=0 'initialize a series to store the occurrences of rejecting the wrong hypothesis

pcurveb2h.xy 'plot the power of the t test curve column 2 (c2) against the test values in column 1 (c1)
!mintv=0 'minimum test value
!maxtv=6 'maximum test value
!ntestv=101 'number of test values
!step=(!maxtv-!mintv)/!ntestv 'step size
'create a !ntestv by 2 matrix called ptest
'the first column stores the wrong test values while the second stores the power of tests
matrix(!ntestv,2) pcurveb3h
series rejecth=0 'initialize a series to store the occurrences of rejecting the wrong hypothesis

for !i=1 to !ntestv
pcurveb3(!i,1)=!mintv+(!i-1)*!step 'the test value
reject=@recode(((b3hat-pcurveb3(!i,1))/seb3hat<tlb) or ((b3hat-pcurveb3(!i,1))/seb3hat>tub), 1, 0 )
pcurveb3(!i,2)=@sum(reject)/!s
pcurveb3wls(!i,1)=!mintv+(!i-1)*!step 'the test value
rejectwls=@recode(((b3hatwls-pcurveb3wls(!i,1))/seb3hatwls<tlb) or ((b3hatwls-pcurveb3wls(!i,1))/seb3hatwls>tub), 1, 0 )
pcurveb3wls(!i,2)=@sum(rejectwls)/!s
pcurveb3h(!i,1)=!mintv+(!i-1)*!step 'the test value
rejecth=@recode(((b3hath-pcurveb3h(!i,1))/seb3hath<tlb) or ((b3hath-pcurveb3h(!i,1))/seb3hath>tub), 1, 0 )
pcurveb3h(!i,2)=@sum(rejecth)/!s
next
'this section plots power curves
pcurveb3.xy 'plot the power of the t test curve column 2 (c2) against the test values in column 1 (c1)
pcurveb3wls.xy 'plot the power of the t test curve column 2 (c2) against the test values in column 1 (c1)
pcurveb3h.xy 'plot the power of the t test curve column 2 (c2) against the test values in column 1 (c1)