# Monte Carlo Applications

- Small Sample Distribution of ADF test (simadf.prg)

## Small sample distribution of ADF tests

The program simadf.prg follows the analysis in Rudebusch (1993) "The Uncertain Unit Root in Real GNP," American Economic Review, 83, 1, 264-272. The first part of the program fits a trend stationary and difference stationary model to log US GDP. The impulse response functions from the two models are plotted and contrasted.

The second part of the program simulates the distribution of the t-statistic for the unit root test under both models. The results are first stored in a matrix object. Then the results are exported to a second workfile where the matrix results are converted to series for further processing. In particular, we plot the kernel density estimates of the simulated distributions of the t-statistic and compute the size and power of the ADF test.

    ' program to simulate the small sample distribution of ADF test
    ' for version 4
    ' last checked 11/2/2000 h
    '
    ' Rudebusch (1993) "The Uncertain Unit Root in Real GNP," AER,
    ' 83, 1, 264-272.

    ' set up workfile
    workfile adftest q 1948:1 1988:4

    ' fill GDP series
    series gdp
    gdp.fill 1284,1295.69995117,1303.80004883,1316.40002441,1305.30004883,1302,1312.59997559,1301.90002441,1350.90002441,1393.5,1445.19995117,1484.5,1504.09997559,1548.30004883,1585.40002441,1596,1607.69995117,1612.09997559,1621.90002441,1657.80004883,1687.30004883,1695.30004883,1687.90002441,1671.19995117,1660.80004883,1658.40002441,1677.69995117,1698.30004883,1742.5,1758.59997559,1778.19995117,1793.90002441,1787,1798.5,1802.19995117,1826.59997559,1836.40002441,1834.80004883,1851.19995117,1830.5,1790.09997559,1804.40002441,1840.90002441,1880.90002441,1904.90002441,1937.5,1930.80004883,1941.90002441,1976.90002441,1971.69995117,1973.69995117,1961.09997559,1977.40002441,2006,2035.19995117,2076.5,2103.80004883,2125.69995117,2142.60009766,2140.19995117,2170.89990234,2199.5,2237.60009766,2254.5,2311.10009766,2329.89990234,2357.39990234,2364,2410.10009766,2442.80004883,2485.5,2543.80004883,2596.80004883,2601.39990234,2626.10009766,2640.5,2657.19995117,2669,2699.5,2715.10009766,2752.10009766,2796.89990234,2816.80004883,2821.69995117,2864.60009766,2867.80004883,2884.5,2875.10009766,2867.80004883,2859.5,2895,2873.30004883,2939.89990234,2944.19995117,2962.30004883,2977.30004883,3037.30004883,3089.69995117,3125.80004883,3175.5,3253.30004883,3267.60009766,3264.30004883,3289.10009766,3259.39990234,3267.60009766,3239.10009766,3226.39990234,3154,3190.39990234,3249.89990234,3292.5,3356.69995117,3369.19995117,3381,3416.30004883,3466.39990234,3525,3574.39990234,3567.19995117,3591.80004883,3707,3735.60009766,3779.60009766,3780.80004883,3784.30004883,3807.5,3814.60009766,3830.80004883,3732.60009766,3733.5,3808.5,3860.5,3844.39990234,3864.5,3803.10009766,3756.10009766,3771.10009766,3754.39990234,3759.60009766,3783.5,3886.5,3944.39990234,4012.10009766,4089.5,4144,4166.39990234,4194.20019531,4221.79980469,4254.79980469,4309,4333.5,4390.5,4387.70019531,4412.60009766,4427.10009766,4460,4515.29980469,4559.29980469,4625.5,4655.29980469,4704.79980469,4734.5,4779.70019531

    ' work with log transformation
    series ly = log(gdp)

    ' estimate trend stationary model
    equation eq_ts.ls ly c @trend+1 ly(-1) ly(-2)

    ' estimate difference stationary model
    equation eq_ds.ls d(ly) c d(ly(-1))

    ' compute impulse response from each model
    ' initialize impulse and response series
    series u = 0
    series ly_ts = 0
    series ly_ds = 0

    ' create the impulse series U
    smpl 48:3 48:3
    series u = 1

    smpl 48:3 88:4

    ' compute impulse response from trend stationary model
    ly_ts = eq_ts.@coefs(3)*ly_ts(-1) + eq_ts.@coefs(4)*ly_ts(-2) + u

    ' compute impulse response from difference stationary model
    ly_ds = ly_ds(-1) + eq_ds.@coefs(2)*d(ly_ds(-1)) + u

    ' plot the impulse responses
    smpl 48:3 58:3
    graph gh1.line ly_ts ly_ds
    ' label graph
    gh1.setelem(1) legend(Trend stationary model)
    gh1.setelem(2) legend(Difference stationary model)
    gh1.addtext(t) Impulse Response Functions
    show gh1

    '--------------------------------------------------------------------------
    ' simulate ADF t-stats from each model
    '--------------------------------------------------------------------------

    ' initialize simulated series
    smpl 48:1 48:2
    series lyts = ly
    series lyds = ly

    ' declare matrix to store simulated ADF t-stats
    !n = 1000				' number of replications
    matrix(!n,2) adft

    ' declare working equation
    equation eq_test

    ' monte carlo loop
    rndseed 123456789
    smpl 48:3 @last
    for !i = 1 to !n
    	' trend stationary model
    	' simulate series
    	lyts = eq_ts.@coefs(1) + eq_ts.@coefs(2)*(@trend+1) + eq_ts.@coefs(3)*lyts(-1) + eq_ts.@coefs(4)*lyts(-2) + eq_ts.@se*nrnd
    	' run ADF test
    	eq_test.ls d(lyts) lyts(-1) c @trend+1 d(lyts(-1))
    	' store t-stat
    	adft(!i,1) = eq_test.@tstat(1)

    	' difference stationary model
    	' simulate series
    	lyds = lyds(-1) + eq_ds.@coefs(1) + eq_ds.@coefs(2)*d(lyds(-1)) + eq_ds.@se*nrnd
    	' run ADF test
    	eq_test.ls d(lyds) lyds(-1) c @trend+1 d(lyds(-1))
    	' store t-stat
    	adft(!i,2) = eq_test.@tstat(1)
    next

    ' store actual ADF t-stat
    eq_test.ls d(ly) ly(-1) c @trend+1 d(ly(-1))
    !tadf = eq_test.@tstat(1)

    ' store results to database
    db db_mcarlo
    store adft

    ' close workfile
    close adftest

    '--------------------------------------------------------------------------
    ' process results as series in new workfile
    '--------------------------------------------------------------------------

    ' create new workfile
    workfile result u 1 !n

    ' fetch results
    fetch adft

    ' declare series
    series t_ts
    series t_ds

    ' convert simulated t-stats to series
    group gr1 t_ts t_ds
    mtos(adft,gr1)

    ' determine knots for kernel density estimates
    ' to ensure it fits in current workfile
    !knots = 100
    if (!knots>!n) then
    	!knots = !n
    endif

    smpl 1 !knots
    for %mod ts ds
    	' store kernel density estimates
    	do t_{%mod}.kdensity(!knots, o=kdmat_{%mod})
    	' convert density estimates to series
    	series kd{%mod}_x
    	series kd{%mod}_f
    	group g{%mod} kd{%mod}_x kd{%mod}_f
    	mtos(kdmat_{%mod}, g{%mod})
    next

    'compute empirical size of test
    ' area of H0 to the left of !tadf
    smpl @all if t_ds<=!tadf
    !siz = @obs(t_ds)/!n

    'compute empirical power of test
    ' area of H1 to the left of !tadf
    smpl @all if t_ts<=!tadf
    !pow = @obs(t_ts)/!n

    ' plot the two density estimates in one graph
    smpl 1 !knots

    group g3 gts gds
    freeze(graph1) g3.xyline(b)
    'graph1.setelem(1) legend(Alternative (trend stationary))
    'graph1.setelem(2) legend(Null (difference stationary))
    graph1.addtext(t) Kernel density estimates of simulated t-distribution
    graph1.legend -display
    graph1.addtext(1.2,0.1) H1 (trend stationary)
    graph1.addtext(1.8,0.8) H0 (difference stationary)
    ' draw vertical dashline at actual t-stat
    graph1.draw(dashline, bottom, rgb(156,156,156)) !tadf
    graph1.addtext(2.8,2.4) size = !siz
    graph1.addtext(2.8,2.6) power = !pow
    show graph1