Method of Moments approach

[BinhPham]

Dear Mr Tom,

Could I have a sample showing how to declare moment conditions for GMM estimation?

I am confusing with the UG hence if you have an example it is much better. The idea is that I want to estimate indivisible RBC's params as in Craig Burnside's matlab codes.

Thanks in advance.

[TomDoan]

What are the moment conditions written as mathematical formulas?

The moment conditions are typically either the Kroneker product of some non-linear u times instruments Z (so you would do FRML's for the u's and an INSTRUMENTS instruction for the Z's) or just direct translations of moment conditions (FRML's for those and CONSTANT as the only "instrument"). An example of the first is GRN7P333.RPF (Greene textbook examples) and of the latter is HAMP410.RPF (Hamilton textbook examples)

[BinhPham]

Dear Tom,

Thank for your guide. I have an issue when specifying moment conditions.
Taking indivisible RBC as an example, it has 13 moment conditions as in Craig Burnside (1999) notes:

Mom1: E[theta - alpha*Y_t/(N_t*C_t)] = 0
Mom2: E[log(A_t) - log(A)*(1-rhoA) - rhoA*log(A_{t-1})] = 0
Mom3: E[[log(A_t) - log(A)*(1-rhoA) - rhoA*log(A_{t-1})]*log(A_{t-1})] = 0
Mom4: E[log(A_t) - log(A)*(1-rhoA) - rhoA*log(A_{t-1}) - sigmaA^2] = 0
Mom5: E[log(Y_t) - Ay - lnGammaX)*t] = 0
Mom6: E[[log(Y_t) - Ay - lnGammaX)*t]*t/T] = 0
Mom7: E[delta - 1 + (K_{t+1} - I_t)/K_t] = 0
Mom8: E[1 - beta*(C_t/C_{t+1})*[(1-alpha)Y_{t+1}/K_{t+1} + (1-delta)]
Mom9: E[hpy_t^2 - sigma_y^2]
Mom10: E[hpc_t^2 - (sigma_c/sigma_y)^2*hpy_t^2] = 0
Mom11: E[hpi_t^2 - (sigma_i/sigma_y)^2*hpy_t^2] = 0
Mom12: E[hpn_t^2 - (sigma_n/sigma_y)^2*hpy_t^2] = 0
Mom13: E[hpn_t^2 - (sigma_n/sigma_yn)^2*hpyn_t^2] = 0

where parameters are: Theta, lnA, rhoA, sigmaA, Ay, lnGammX, Delta, Alpha, sigma_y, and four ratio sigma_.../sigma_y.

Observables are: C I Y N K and corresponding hp filtered of log(.): hpc hpi hpy hpn hpk, and hpyn = hpy - hpn
Writing moment in frml is trivial as I can write:

Code: Select all

nonlin Theta lnA rhoA sigmaA Ay lnGammX Delta Alpha sigmaY sigmaCY sigmaIY sigmaNY sigmaNYN
* setting initial paramater values
compute Theta = ...
...
frml mom1 = theta - alpha*Y/(C*N)
frml mom2 = logA - logA*(1-rhoA) - rhoA*logA{1}
...

But Solow residual A_t is not observed, hence, it is estimated before the NLSYSTEM command as:

Code: Select all

set trend = t
set logA   = log(Y) - log(K)*(1-alpha) - log(N)*alpha - trend*lnGammX*alpha

The problem is that series logA depends upon lnGammX and alpha which are the parameters would be estimated.
Since we need to set up initial values of 13 parameters so the series logA will be computed at the outset but it will not be updated during the optimization. In the Burnside (1999) matlab code, the function gamerr.m has a line as follows (with my annotations):

Code: Select all

function mom=gmmerr(b)
   global xdata;
   loga = log(xdata(1:115,var_y))-log(xdata(1:115,var_k))*(1-b(8,1)) ...
             - log(xdata(1:115,var_n))*b(8,1)-t0*b(6,1)*b(8,1) ;
   ...
   %Moment 2: E[log(A_t) - log(A)*(1-rhoA) - rhoA*log(A_{t-1})] = 0
   mom1 = ...
   mom2 = loga(2:114,1)-b(2,1)-loga(1:113,1)*b(3,1) ;
   ...
   mom13= ...
mom = [mom1 mom2 ... mom13]; 

How can I set up the dynamic series logA for GMM estimation? I think it is the need of a helper function. Could you please show me a tip?

Thank you and best regards,

[TomDoan]

Define a FRML which gives you log A(t) from the other stuff:

frml logAF = log(Y) - log(K)*(1-alpha) - log(N)*alpha - trend*lnGammX*alpha

then your mom2 (for instance) becomes

frml mom2 = logAF{0}-log(A)*(1-rhoA)-rhoA*logAF{1}

[BinhPham]

Dear Tom,

Thank for your help. I code it up but the estimation seems to be not convergence even I set nlsystem(...,iterations=10000), and the estimates are bad. Another question is should I include frml logA = ... in nlsystem ? I estimate both case but the results are far from Burnside (1999) matlab code.

Could you please give me a hint?

The target results are:
Theta 5.1199
lnA 0.1936
rhoA 0.9663
sigA 0.0119
Ay 8.5722
lnGamm 0.0041
Delta 0.0209
Alpha 0.6553
SigY 0.0192
SigCY 0.4366
SigIY 2.2240
SigNY 0.8585
SigNYN 1.2210

My bad results are:

Code: Select all

GMM-Factored Weight Matrix
NO CONVERGENCE IN 500 ITERATIONS
LAST CRITERION WAS  0.0001080
TRY INCREASING ITERS OPTION

Quarterly Data From 1955:03 To 1984:01
Usable Observations                          113
Skipped/Missing (from 115)                     2
Function Value                      112.76075070
J-Specification(1)                      112.7608
Significance Level of J                0.0000000

    Variable                          Coeff       Std Error      T-Stat      Signif
***************************************************************************************
1.  THETA                               4.513095     0.051122     88.28005  0.00000000
2.  LNA                              -317.529353  1382.173148     -0.22973  0.81830005
3.  RHOA                                0.999594     0.001741    574.26025  0.00000000
4.  SIGA                               -0.011397     0.011744     -0.97046  0.33181815
5.  AY                                 -0.079760     0.048126     -1.65731  0.09745694
6.  LNGAMM                              0.105594     0.003131     33.72166  0.00000000
7.  DELTA                               0.021604     0.000096    225.84222  0.00000000
8.  ALPHA                               0.642989     0.004735    135.79417  0.00000000
9.  SIGY                                0.012660     0.000843     15.01289  0.00000000
10. SIGCY                               0.504926     0.036085     13.99260  0.00000000
11. SIGIY                               2.058394     0.086683     23.74612  0.00000000
12. SIGNY                               0.926871     0.091827     10.09371  0.00000000
13. SIGNYN                              1.001755     0.073850     13.56475  0.00000000

My code is

Code: Select all

* Housekeeping work, reset the current environment
END(RESET)
*
* Loading dataset of Burnside (1999) note, i.e. hoarding.dat dataset of Christiano and Eichenbaum (CE) (1992).
* Variables from left to right as in the paper are per-capita levels of C I G Y Hh K and He, respectively.
* Note that: Hh and He are the household and establishment hours series as in CE (1992).
* Data need to be filtered (HP) before estimation.
OPEN DATA "hoarding.dat"
CALENDAR(Q) 1955:03
DATA(FORMAT=PRN,NOLABELS,ORG=COLUMNS) 1955:03 1984:01 C DK G Y Hh K He

* Looing at descriptive statistics
TABLE
set trend = t
set n  = Hh/1369
set lc = log(C)
set li = log(DK)
set ly = log(Y)
set lk = log(K)
set ln = log(N)
compute lambda = 100
*
*The sample has 115 data points but minus the first and the last obs
*since endogenous vars involve lag(1) and lead(1)
*
compute NOBS   = %allocend()-2

* HP filter data to get cyclical components
filter(type=hp,tuning=lambda) lc / lchpt
filter(type=hp,tuning=lambda) li / lihpt
filter(type=hp,tuning=lambda) ly / lyhpt
filter(type=hp,tuning=lambda) ln / lnhpt
set hpc  = lc - lchpt
set hpi  = li - lihpt 
set hpy  = ly - lyhpt 
set hpn  = ln - lnhpt
set hpyn = hpy - hpn 

SPGRAPH(hfields=2,vfields=2,window="HP filtered data")
   GRAPH(STYLE=LINE,HEADER="C")
   # hpc /
   GRAPH(STYLE=LINE,HEADER="I")
   # hpi /
   GRAPH(STYLE=LINE,HEADER="Y")
   # hpy /
   GRAPH(STYLE=LINE,HEADER="N")
   # hpn /
SPGRAPH(done)

*Parameters: Theta,lnA,rhoA,sigA,Ay,lnGamm,Delta,Alpha, and five Data Moments
*
nonlin theta lnA rhoA sigA Ay lnGamm Delta Alpha sigY sigCY sigIY sigNY sigNYN
*
*Setting initial values
*bstart=[ 1.0 1.0 0.9 0.01 1.0 0.01 0.02 0.6 0.01 1.0 1.0 1.0 1.0 ]' ;
*
compute beta_  = 1.03^(-0.25)
compute Theta  = 1.0
compute lnA    = 1.0
compute rhoA   = 0.9
compute sigA   = 0.01
compute Ay     = 1.0
compute lnGamm = 0.01
compute Delta  = 0.02
compute Alpha  = 0.6
compute sigY   = 0.01
compute sigCY  = 1.0
compute sigIY  = 1.0
compute sigNY  = 1.0
compute sigNYN = 1.0

frml logA      = log(Y) - log(K)*(1-alpha) - log(N)*alpha -trend*lnGamm*alpha
frml mom1      = theta - alpha*Y/(C*N)
frml mom2      = 1 - (K{-1} - DK)/K - delta
frml mom3      = 1 - beta_*(C/C{-1})*((1-alpha)*Y{-1}/K{-1} + (1-delta))
frml mom4      = log(Y) - Ay - trend*lnGamm
frml mom5      = (log(Y) - Ay - trend*lnGamm)*trend/NOBS
frml mom6      = logA{0} - lnA*(1-rhoA) - rhoA*logA{1}
frml mom7      = (logA{0} - lnA*(1-rhoA) - rhoA*logA{1})*logA{1}
frml mom8      = (logA{0} - lnA*(1-rhoA) - rhoA*logA{1})^2 - sigA^2
frml mom9      = hpy^2 - sigY^2
frml mom10     = hpc^2 - sigCY^2*hpy^2
frml mom11     = hpi^2 - sigIY^2*hpy^2
frml mom12     = hpn^2 - sigNY^2*hpy^2
frml mom13     = hpn^2 - sigNYN^2*hpyn^2

instruments constant
nlsystem(instruments, iterations = 500) / logA mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 mom9 mom10 mom11 mom12 mom13

[TomDoan]

No. The log A FRML is a definition, not a moment condition. If you feed in the Burnside's estimates, your model agrees, so the problem is bad guess values. Offhand, it looks like your Ay (at minimum) is way off. You should be able to guess AY and lnGamm by just running a linear regression of log Y on constant and trend. (Since it's a just identified model, those "guess" values have to be the final estimates).

[BinhPham]

Dear Tom,

Thank for your suggestion. At the end of the day, I only need the better initial values.
The results now sound good with lnA=0.2, rhoA=0.95 and Ay, LnGamm are set to OLS log(Y) over const trend.
There is only lnA=0.2 can give the reproducible result though my estimate lnA is bigger than Burnside. I do know why setting higher lnA will cause ridiculous outputs. Perhaps, internal optimization matters.
Best,

Code: Select all

GMM-Factored Weight Matrix
Convergence in   169 Iterations. Final criterion was  0.0000022 <=  0.0000100

Quarterly Data From 1955:03 To 1984:01
Usable Observations                          113
Skipped/Missing (from 115)                     2
Function Value                        0.00000180

    Variable                          Coeff       Std Error      T-Stat      Signif
***************************************************************************************
1.  THETA                           5.1198536604 0.0412155315    124.22146  0.00000000
2.  LNA                             5.7360535670 0.0687946008     83.37941  0.00000000
3.  RHOA                            0.9662726368 0.0262999796     36.74043  0.00000000
4.  SIGA                            0.0118984342 0.0008937807     13.31248  0.00000000
5.  AY                              8.5680945907 0.0094475481    906.91198  0.00000000
6.  LNGAMM                          0.0040852240 0.0001553471     26.29739  0.00000000
7.  DELTA                           0.0208688241 0.0001347813    154.83477  0.00000000
8.  ALPHA                           0.6552578566 0.0047914672    136.75516  0.00000000
9.  SIGY                            0.0122149612 0.0008767491     13.93211  0.00000000
10. SIGCY                           0.4189475114 0.0335720930     12.47904  0.00000000
11. SIGIY                           2.0654256739 0.0870509453     23.72663  0.00000000
12. SIGNY                           0.9702197280 0.0913483730     10.62109  0.00000000
13. SIGNYN                          1.0220311981 0.0737498016     13.85809  0.00000000

This is the final code for whom wants to do GMM estimation of indivisible RBC model. This code is based on Burnside (1999) matlab codes. Data file is hoarding.dat as above.

Code: Select all

* Housekeeping work, reset the current environment
END(RESET)
*
* Loading dataset of Burnside (1999) note, i.e. hoarding.dat dataset of Christiano and Eichenbaum (CE) (1992).
* Variables from left to right as in the paper are per-capita levels of C I G Y Hh K and He, respectively.
* Note that: Hh and He are the household and establishment hours series as in CE (1992).
* Data need to be filtered (HP) before estimation.
OPEN DATA "hoarding.dat"
CALENDAR(Q) 1955:03
DATA(FORMAT=PRN,NOLABELS,ORG=COLUMNS) 1955:03 1984:01 C DK G Y Hh K He

* Looing at descriptive statistics
TABLE
set trend = t
set n  = Hh/1369
set lc = log(C)
set li = log(DK)
set ly = log(Y)
set lk = log(K)
set ln = log(N)
compute lambda = 100
*
*The sample has 115 data points but minus the first and the last obs
*since endogenous vars involve lag(1) and lead(1)
*
compute NOBS   = %allocend()-2

* HP filter data to get cyclical components
filter(type=hp,tuning=lambda) lc / lchpt
filter(type=hp,tuning=lambda) li / lihpt
filter(type=hp,tuning=lambda) ly / lyhpt
filter(type=hp,tuning=lambda) ln / lnhpt
set hpc  = lc - lchpt
set hpi  = li - lihpt
set hpy  = ly - lyhpt
set hpn  = ln - lnhpt
set hpyn = hpy - hpn

SPGRAPH(hfields=2,vfields=2,window="Observables")
   GRAPH(STYLE=LINE,HEADER="Y")
   # Y /
   GRAPH(STYLE=LINE,HEADER="C")
   # C /
   GRAPH(STYLE=LINE,HEADER="I (DK)")
   # DK /
   GRAPH(STYLE=LINE,HEADER="N")
   # N /
SPGRAPH(done)

SPGRAPH(hfields=2,vfields=2,window="HP filtered data")
   GRAPH(STYLE=LINE,HEADER="C")
   # hpc /
   GRAPH(STYLE=LINE,HEADER="I")
   # hpi /
   GRAPH(STYLE=LINE,HEADER="Y")
   # hpy /
   GRAPH(STYLE=LINE,HEADER="N")
   # hpn /
SPGRAPH(done)
*
*Parameters: Theta,lnA,rhoA,sigA,Ay,lnGamm,Delta,Alpha, and five Data Moments
*
nonlin theta lnA rhoA sigA Ay LnGamm Delta Alpha sigY sigCY sigIY sigNY sigNYN
*
*Setting initial values
*bstart=[ 1.0 1.0 0.9 0.01 1.0 0.01 0.02 0.6 0.01 1.0 1.0 1.0 1.0 ]' ;
*
* Auxiliary OLS to get Ay and LnGamm 
*
linreg ly 1955:4 1983:4
# constant trend
compute Ay     = %beta(1)
compute lnGamm = %beta(2)

compute beta_  = 1.03^(-0.25)
compute Theta  = 1
compute lnA    = 0.2
compute rhoA   = 0.95
compute sigA   = 0.01
compute Delta  = 0.02
compute Alpha  = 0.6
compute sigY   = 0.01
compute sigCY  = 1.0
compute sigIY  = 1.0
compute sigNY  = 1.0
compute sigNYN = 1.0

frml logA      = log(Y) - log(K)*(1-alpha) - log(N)*alpha -trend*lnGamm*alpha
frml mom1      = theta - alpha*Y/(C*N)
frml mom2      = 1 - (K{-1} - DK)/K - delta
frml mom3      = 1 - beta_*(C/C{-1})*((1-alpha)*Y{-1}/K{-1} + (1-delta))
frml mom4      = log(Y) - Ay - trend*lnGamm
frml mom5      = (log(Y) - Ay - trend*lnGamm)*trend/NOBS
frml mom6      = logA{0} - lnA*(1-rhoA) - rhoA*logA{1}
frml mom7      = (logA{0} - lnA*(1-rhoA) - rhoA*logA{1})*logA{1}
frml mom8      = (logA{0} - lnA*(1-rhoA) - rhoA*logA{1})^2 - sigA^2
frml mom9      = hpy^2 - sigY^2
frml mom10     = hpc^2 - sigCY^2*hpy^2
frml mom11     = hpi^2 - sigIY^2*hpy^2
frml mom12     = hpn^2 - sigNY^2*hpy^2
frml mom13     = hpn^2 - sigNYN^2*hpyn^2

instruments constant
nlsystem(instruments, iterations = 500) / mom1 mom2 mom3 mom4 mom5 mom6 mom7 mom8 mom9 mom10 mom11 mom12 mom13

[TomDoan]

The problem is that the log A is borderline non-stationary depending upon the values of alpha and lnGamm. The particular way that the log A equation is formulated (with lnA as the mean) can cause some numerical problems if rhoA ventures above 1. This is also an exactly identified model, which means that the solution is the same regardless of the weight matrix, but different methods of doing the weight matrix can affect how easily it finds the solution. Possibly just putting in SWMATRIX=%IDENTITY(13) will let it solve more easily from a wide variety of guess values. However, there's really no excuse in a situation like this from not computing reasonable guesses to start---only the log A isn't observable.