LM test for mean stationarity
How to implement the LM test for mean stationarity in "Testing Initial Conditions in Dynamic Panel Data Models" with G. Calzolari (Econometric Reviews, 2020)
The command xtabond2 in STATA allows recovering the matrices that are needed for computation of the LM test suggested by Magazzini & Calzolari when the option svmat is specified. Also you have to use the Mata implementation in Mata's speed-favoring mode:
. mata: mata set matafavor speedIn the pure dynamic case:
. xtabond2 y L.y, gmmstyle(L.y) twostep robust svmat noconst . mat bbb2=e(b)You will need to compute the matrix G, and add the column that contains the first derivative of the moment condition with respect to the additional parameter . Set a local variable TT equal to the number of time periods.
These will be equal to zero for the moment conditions related to the first-differenced equations, and equal to "-bsys2^(`tt'-2-1)" for the moment conditions related to the level equations.
Denote TT the maximum number of time periods.
After the command xtabond2, you can type:
. mat G=-(e(Z))'*(e(X)) /* first derivative of the moment conditions with respect to beta */. mat gpsi =J(colsof(e(Z)),1,0) /* the additional column related to the additional parameter */. mat bsys2=e(b). scalar bsys2=bsys2[1,1]. forvalues tt=3/`TT' { mat gpsi[colnumb(e(Z), "Levels eq:LD.y/"+strofreal(`tt')),1]=-bsys2^(`tt'-2-1) }. mat G=(G,gpsi). mat testcm = e(Ze)'*e(A2)*G*invsym(G'*e(A2)*G)*G'*e(A2)*e(Ze)When an endogenous regressor x is included, the number of columns to be added to G is 1 (related to the parameter \psi) plus T-1 (related to \xi_1,...,\xi_(T-1)):
. xtabond2 y L.y x, gmmstyle(y x, lag(2 .)) twostep robust svmat noconst h(2) mat bsys2=e(b). scalar bly=bsys2[1,1]. scalar bx=bsys2[1,2]. mat G=-(e(Z))'*(e(X)). * parameter psi related to E(\Delta y_{it-1} \epsilon_{it}). mat gpsi =J(colsof(e(Z)),1,0). forvalues ii=3/`TT' { mat gpsi[colnumb(e(Z), "Levels eq:LD.y/"+strofreal(`ii')),1]=-bly^(`ii'-2-1) }. * parameters xi(1)...xi(T-1) related to E(\Delta x_{i1} \epsilon_{i2}),...,E(\Delta x_{i,T-1} \epsilon_{iT}). mat gxi = J(colsof(e(Z)),`TT'-1-1,0). forvalues ii=3/`TT' { mat gxi[colnumb(e(Z), "Levels eq:LD.x/"+strofreal(`ii')),(`ii'-2)]=-1 }. forvalues ii=4/`TT' { local ii4=`ii'-4 forvalues ss=0/`ii4' { local jj = `ii'-`ss' mat gxi[colnumb(e(Z), "Levels eq:LD.y/"+strofreal(`ii')),(`jj'-2)]=-bx*bly^(`ss') } }. mat G=(G,gpsi,gxi). mat testcm = e(Ze)'*e(A2)*G*invsym(G'*e(A2)*G)*G'*e(A2)*e(Ze)
Example: Data from Ziliak (1997)
We will use data from Ziliak (1997) in order to understand how to implement the proposed test procedure on real data. The dataset is also used in the book by Cameron & Trivedi (2005), and the dataset can be downloaded at the following website.
We will consider a pure dynamic model with no regressors, as well as the cases of strict exogenous regressor, predetermined x, and endogenous variables. In all cases, before applying the proposed testing procedure, you should check whether the assumption for the application of the difference-GMM estimator are satisfied. The LM test procedure we propose only checks the validity of the mean stationarity assumption, needed for consistency of the system-GMM estimator.
Dynamic model with no regressor
First consider the model:
lnhr[i,t] = beta0 + beta1 lnhr[i,t-1] + tau[t] + u[i] + e[i,t] The system-GMM estimator can be obtained by using:
. xtabond2 lnhr l.lnhr dyear3-dyear10, gmmstyle(l.lnhr) ivstyle(dyear3-dyear10, eq(level)) twostep robust svmatThe option "svmat" ensures that all matrices that are needed for computation are stored into memory.
In order to compute the test, we need to "augment" the moment conditions related to the equation in levels. Before doing that we store the estimated coefficient in the vector "bsys2" and we compute the matrix G:
. mat bsys2=e(b). mat G=-(e(Z))'*(e(X))The matrix G contains the first derivative of the moment condition with respect to the parameter of interest. In order to compute the LM test, we have to "augment" the moment conditions for the level equation, and compute the first derivative with respect to this additional parameter. The vector of first derivative will be then included in the matrix G.
. mat gpsi =J(colsof(e(Z)),1,0). mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1981"),1]=-bsys2[1,1]^0. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1982"),1]=-bsys2[1,1]^1. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1983"),1]=-bsys2[1,1]^2. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),1]=-bsys2[1,1]^3. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),1]=-bsys2[1,1]^4. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),1]=-bsys2[1,1]^5. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),1]=-bsys2[1,1]^6. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),1]=-bsys2[1,1]^7. mat G=(G,gpsi)The LM test can be now computed as:
. mat testcm = e(Ze)'*e(A2)*G*invsym(G'*e(A2)*G)*G'*e(A2)*e(Ze)The test is chi-squared distributed with 1 degree of freedom (equal to the number of additional parameters).
The case of strictly exogenous regressors
We will now include the variable "lnwg" in the equation and treat it as a strictly exogenous variable. The system-GMM equation can be obtained by:
. xtabond2 lnhr l.lnhr lnwg dyear3-dyear10, gmmstyle(l.lnhr) ivstyle(dyear3-dyear10, eq(level)) ivstyle(lnwg, eq(diff)) robust twostep svmatIn this case, there are no additional moment conditions to be tested that are related to the strictly exogenous variable. The computation of the test can proceed in the same way as in the pure dynamic case.
Treating lnwg as predetermined variable
If we want to treat lnwg as predetermined variable, then we need to specify the variable in the gmmstyle option:
. xtabond2 lnhr l.lnhr lnwg dyear3-dyear10, gmmstyle(l.lnhr lnwg) ivstyle(dyear3-dyear10, eq(level)) robust twostep svmat. mat bsys2=e(b). mat G=-(e(Z))'*(e(X))As for the moment conditions related to the level of the dependent variable, the construction of the gpsi matrix is unaffected:
. mat gpsi =J(colsof(e(Z)),1,0). mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1981"),1]=-bsys2[1,1]^0. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1982"),1]=-bsys2[1,1]^1. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1983"),1]=-bsys2[1,1]^2. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),1]=-bsys2[1,1]^3. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),1]=-bsys2[1,1]^4. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),1]=-bsys2[1,1]^5. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),1]=-bsys2[1,1]^6. mat gpsi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),1]=-bsys2[1,1]^7However, we will now need to "augment" also the moment conditions related to lnwg. As the variable is treated as predetermined, we have 9 additional moment conditions: E(D.lnwg(80) eps(80)) = 0 ,..., E(D.lnwg(88) eps(88)) = 0 (with D the first differencing operator). We will need to "augment" all these moment conditions by letting E(D.lnwg(80) eps(80)) - PW1 = 0 ,..., E(D.lnwg(88) eps(88)) -PW9 = 0 (with PW1,...,PW9 representing the additional parameters).
In this case, the mean stationarity assumption will be satisfied when PSI = PW1 = ... PW9 = 0. We will need to compute the first derivative also with repsect to those additional parameters:
. mat gxi=J(colsof(e(Z)),9,0). mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1980"),1]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1981"),2]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1982"),3]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1983"),4]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1984"),5]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1985"),6]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1986"),7]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1987"),8]=-1. mat gxi[colnumb(e(Z), "Levels eq:D.lnwg/1988"),9]=-1These parameters also enter the level moment conditions related to y:
. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1982"),3]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1983"),3]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),3]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),3]=-bsys2[1,1]^3*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),3]=-bsys2[1,1]^4*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),3]=-bsys2[1,1]^5*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),3]=-bsys2[1,1]^6*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1983"),4]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),4]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),4]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),4]=-bsys2[1,1]^3*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),4]=-bsys2[1,1]^4*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),4]=-bsys2[1,1]^5*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),5]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),5]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),5]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),5]=-bsys2[1,1]^3*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),5]=-bsys2[1,1]^4*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),6]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),6]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),6]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),6]=-bsys2[1,1]^3*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),7]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),7]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),7]=-bsys2[1,1]^2*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),8]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),8]=-bsys2[1,1]*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),9]=-bsys2[1,2] . mat G=(G,gpsi,gxi). mat testcm = e(Ze)'*e(A2)*G*invsym(G'*e(A2)*G)*G'*e(A2)*e(Ze)The test is chi-squared distributed with 10 degrees of freedom.
Treat lnwg as an endogenous variables
In the case the variable lnwg is treated as an endogenous variable, lnwg(t-1) is no longer a valid instrument: The system-GMM estimator can be obtained as:
. xtabond2 lnhr l.lnhr lnwg dyear3-dyear10, gmmstyle(l.lnhr) gmmstyle(lnwg, lag(2 .)) ivstyle(dyear3-dyear10, eq(level)) robust twostep svmat. mat bsys2=e(b). mat G=-(e(Z))'*(e(X))The matrix gpsi is build as in the previous cases, so we will omit it here.
As for the moment conditions related to lnwg, we now have 8 additional moment conditions, that can be "augmented" by using the parameters PW1,..., PW8:
E(D.lnwg(80) eps(81)) - PW1 =0
E(D.lnwg(80) eps(81)) - PW2 =0
...
E(D.lnwg(87) eps(88)) - PW8 =0
The matrix gxi can therefore be constructed as:
. mat gxi=J(colsof(e(Z)),8,0). mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1981"),1]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1982"),2]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1983"),3]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1984"),4]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1985"),5]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1986"),6]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1987"),7]=-1. mat gxi[colnumb(e(Z), "Levels eq:LD.lnwg/1988"),8]=-1Moreover, the parameter PW1, ..., PW8 also appears in the moment conditions related to y:
. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1982"),2]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1983"),2]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),2]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),2]=-bsys2[1,1]^3*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),2]=-bsys2[1,1]^4*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),2]=-bsys2[1,1]^5*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),2]=-bsys2[1,1]^6*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1983"),3]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),3]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),3]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),3]=-bsys2[1,1]^3*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),3]=-bsys2[1,1]^4*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),3]=-bsys2[1,1]^5*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1984"),4]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),4]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),4]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),4]=-bsys2[1,1]^3*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),4]=-bsys2[1,1]^4*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1985"),5]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),5]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),5]=-bsys2[1,1]^2*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),5]=-bsys2[1,1]^3*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1986"),6]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),6]=-bsys2[1,1]*bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),6]=-bsys2[1,1]^2*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1987"),7]=-bsys2[1,2]. mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),7]=-bsys2[1,1]*bsys2[1,2] . mat gxi[colnumb(e(Z), "Levels eq:LD.lnhr/1988"),8]=-bsys2[1,2]The test is finally computed as:
. mat G=(G,gpsi,gxi). mat testcm = e(Ze)'*e(A2)*G*invsym(G'*e(A2)*G)*G'*e(A2)*e(Ze)Under the null hyopthesis, the testcm is chi-squared distributed with 9 degrees of freedom.
Example: Data from Bond (2002)
We now consider the dataset usbal89.dta used in Bond (2002).
With respect to previous example, when using these data, only lags 3 or older can be used as legitimate instruments, and lagged values of the regressors are included in the equation of interest. The proposed test procedure can still be applied; however, we will need to understand how to "augment" the moment conditions related to the equation in levels.
Open the dataset
. use usbal89and set it as a panel dataset:
. xtset id yearWe build on the models originally studied by Bond (2002). See here for details. However, we will use the twostep option to compute the System GMM estimator
First set mata in "favour-speed" mode -- otherwise some of the matrices that are needed for computation will not be stored in memory.
. mata: mata set matafavor speedThe proposed test procedure can be employed to test the validity of the mean stationarity assumption. Before applying the procedure we need to check whether the difference-GMM estimator is properly specificied:
. xi: xtabond2 y l.y n l.n k l.k i.year , gmm(y n k, lag(3 .)) iv(i.year, equation(level)) twostep robust noleveleqThe diagnostic statistics seem fine. We compute the system-GMM estimator and then compute the LM test for the validity of the mean stationarity assumption:
. xi: xtabond2 y l.y n l.n k l.k i.year , gmm(y n k, lag(3 .)) iv(i.year, equation(level)) twostep robust svmat. mat bsys2=e(b). mat G = -(e(Z))'*(e(X))Before proceeding with the computations, let us explicitly write all the moment conditions in order to understand how to "augment" the parameters in the moment conditions
The equation of interest is
y(i,t) = b1 y(i,t-1) + b2 n(i,t) + b3 n(i,t-1) + b4 k(i,t) + b5 k(i,t-1) + time dummies + eps(i,t)
For the difference-GMM estimator we are considering the following moment conditions:
E(y(i,t-j) D.eps(i,t)) = 0 with j>=3
(where D represent the first-differencing operator)
-- note that we are not using y(i,t-2) as instrument as these do not appear to be legitimate instruments in this application (Bond, 2002).
We will explicitly write all the moment conditions in levels, as these are the ones that need to be "augmented".
We have:
- 5 additional moment conditions (and 5 parameters PN1,...,PN5) related to n
E(D.n(i,83) eps(i,85)) - PN1 = 0
E(D.n(i,84) eps(i,86)) - PN2 = 0
E(D.n(i,85) eps(i,87)) - PN3 = 0
E(D.n(i,86) eps(i,88)) - PN4 = 0
E(D.n(i,87) eps(i,89)) - PN5 = 0
- 5 additional moment conditions (and 5 parameters PK1,...,PK5) related to k
E(D.k(i,83) eps(i,85)) - PK1 = 0
E(D.k(i,84) eps(i,86)) - PK2 = 0
E(D.k(i,85) eps(i,87)) - PK3 = 0
E(D.k(i,86) eps(i,88)) - PK4 = 0
E(D.k(i,87) eps(i,89)) - PK5 = 0
- 5 additional moment conditions (and 1 parameter PSI) related to y
E(D.y(i,83) eps(i,85)) - PSI = 0
E(D.y(i,84) eps(i,86)) - b1*PSI - b2 PN2 - b3 PN1 - b4 PK2 - b5 PK1 = 0
E(D.y(i,85) eps(i,87)) - b1^2*PSI - b1*b2 PN2 - b1*b3 PN1 - b1*b4 PK2 - b1*b5 PK1 - b2 PN3 - b3 PN2 - b4 PK3 - b5 PK2 =
= E(D.y(i,85) eps(i,87)) - b1^2*PSI - b2 PN3 - (b1*b2+b3) PN2 - b1*b3 PN1 - b4 PK3 - (b1*b4+b5) PK2 - b1*b5 PK1 = 0
E(D.y(i,86) eps(i,88)) - b1^3*PSI - b1*b2 PN3 - b1*(b1*b2+b3) PN2 - b1^2*b3 PN1 - b1*b4 PK3 - b1*(b1*b4+b5) PK2 - b1^2*b5 PK1 - b2 PN4 - b3 PN3 - b4 PK4 - b5 PK3 =
= E(D.y(i,86) eps(i,88)) - b1^3*PSI - b2 PN4 - (b1*b2+b3) PN3 - b1*(b1*b2+b3) PN2 - b1^2*b3 PN1 - b4 PK4 - (b1*b4+b5) PK3 - b1*(b1*b4+b5) PK2 - b1^2*b5 PK1 = 0
E(D.y(i,87) eps(i,89)) - b1^4*PSI - b1*b2 PN4 - b1*(b1*b2+b3) PN3 - b1^2*(b1*b2+b3) PN2 - b1^3*b3 PN1 - b1*b4 PK4 - b1*(b1*b4+b5) PK3 - b1^2*(b1*b4+b5) PK2 - b1^3*b5 PK1 - b2 PN5 - b3 PN4 - b4 PK5 - b5 PK4 =
= E(D.y(i,87) eps(i,89)) - b1^4*PSI - b2 PN5 - (b1*b2+b3) PN4 - b1*(b1*b2+b3) PN3 - b1^2*(b1*b2+b3) PN2 - b1^3*b3 PN1 - b4 PK5 - (b1*b4+b5) PK4 - b1*(b1*b4+b5) PK3 - b1^2*(b1*b4+b5) PK2 - b1^3*b5 PK1 = 0
In order to understand how the moment conditions related to y are built, let us consider the formula for D.y(i,84):
D.y(i,84) = b1 D.y(i,83) + b2 D.n(i,84) + b3 D.n(i,83) + b4 D.k(i,84) + b5 D.k(i,83) + D.(time dummies) + D.eps(i,84)
so that,
E[D.y(i,84)*eps(i,86)] =
= b1 E[D.y(i,83)*eps(i,86)] + b2 E[D.n(i,84)*eps(i,86)] + b3 E[D.n(i,83)*eps(i,86)] + b4 E[D.k(i,84)*eps(i,86)] + b5 E[D.k(i,83)*eps(i,86)]
Given the moment conditions in the difference equations, we can write
= b1 E[D.y(i,83)*eps(i,85)] + b2 E[D.n(i,84)*eps(i,86)] + b3 E[D.n(i,83)*eps(i,85)] + b4 E[D.k(i,84)*eps(i,86)] + b5 E[D.k(i,83)*eps(i,85)]
= b1 PSI + b2 PN2 + b3 PN1 + b4 PN2 + b5 PN1
The moment conditions at subsequent time period can be obtained by further substitution.
In order to compute the test, we need to build the first derivative related to the addtional parameters. First let us condier the first derivative related to the parameter PSI:
. mat gpsi =J(colsof(e(Z)),1,0). mat gpsi[colnumb(e(Z), "Levels eq:L2D.y/1985"),1]=-bsys2[1,1]^0. mat gpsi[colnumb(e(Z), "Levels eq:L2D.y/1986"),1]=-bsys2[1,1]^1. mat gpsi[colnumb(e(Z), "Levels eq:L2D.y/1987"),1]=-bsys2[1,1]^2. mat gpsi[colnumb(e(Z), "Levels eq:L2D.y/1988"),1]=-bsys2[1,1]^3. mat gpsi[colnumb(e(Z), "Levels eq:L2D.y/1989"),1]=-bsys2[1,1]^4Now, we compute the first derivative with respect to the parameters PN1,...,PN5 related to moment conditions E(\Delta n_{i,83} \epsilon_{i,84}),...,E(\Delta n_{i,87} \epsilon_{i,89}):
. mat gn = J(colsof(e(Z)),5,0). mat gn[colnumb(e(Z), "Levels eq:L2D.n/1985"),1]=-1. mat gn[colnumb(e(Z), "Levels eq:L2D.n/1986"),2]=-1. mat gn[colnumb(e(Z), "Levels eq:L2D.n/1987"),3]=-1. mat gn[colnumb(e(Z), "Levels eq:L2D.n/1988"),4]=-1. mat gn[colnumb(e(Z), "Levels eq:L2D.n/1989"),5]=-1These parameters also appear in the moment conditions related to y(i,t-1):
. mat gn[colnumb(e(Z), "Levels eq:L2D.y/1986"),1] = -bsys2[1,3]. mat gn[colnumb(e(Z), "Levels eq:L2D.y/1987"),1] = -bsys2[1,1]*bsys2[1,3]. mat gn[colnumb(e(Z), "Levels eq:L2D.y/1988"),1] = -bsys2[1,1]^2*bsys2[1,3]. mat gn[colnumb(e(Z), "Levels eq:L2D.y/1989"),1] = -bsys2[1,1]^3*bsys2[1,3] . mat gn[colnumb(e(Z), "Levels eq:L2D.y/1986"),2]=-bsys2[1,2]. mat gn[colnumb(e(Z), "Levels eq:L2D.y/1987"),2]=-(bsys2[1,1]*bsys2[1,2]+bsys2[1,3]). mat gn[colnumb(e(Z), "Levels eq:L2D.y/1988"),2]=-bsys2[1,1]*(bsys2[1,1]*bsys2[1,2]+bsys2[1,3]). mat gn[colnumb(e(Z), "Levels eq:L2D.y/1989"),2]=-bsys2[1,1]^2*(bsys2[1,1]*bsys2[1,2]+bsys2[1,3]) . mat gn[colnumb(e(Z), "Levels eq:L2D.y/1987"),3]=-bsys2[1,2] . mat gn[colnumb(e(Z), "Levels eq:L2D.y/1988"),3]=-(bsys2[1,1]*bsys2[1,2]+bsys2[1,3]). mat gn[colnumb(e(Z), "Levels eq:L2D.y/1989"),3]=-bsys2[1,1]*(bsys2[1,1]*bsys2[1,2]+bsys2[1,3]) . mat gn[colnumb(e(Z), "Levels eq:L2D.y/1988"),4]=-bsys2[1,2]. mat gn[colnumb(e(Z), "Levels eq:L2D.y/1989"),4]=-(bsys2[1,1]*bsys2[1,2]+bsys2[1,3]) . mat gn[colnumb(e(Z), "Levels eq:L2D.y/1989"),5]=-bsys2[1,2] Then compute the first derivative with respect to the parameters PK1,...,PK5 related to the moment conditions E(\Delta k_{i,83} \epsilon_{i,84}),...,E(\Delta k_{i,87} \epsilon_{i,89}):
. mat gk = J(colsof(e(Z)),5,0). mat gk[colnumb(e(Z), "Levels eq:L2D.k/1985"),1]=-1. mat gk[colnumb(e(Z), "Levels eq:L2D.k/1986"),2]=-1. mat gk[colnumb(e(Z), "Levels eq:L2D.k/1987"),3]=-1. mat gk[colnumb(e(Z), "Levels eq:L2D.k/1988"),4]=-1. mat gk[colnumb(e(Z), "Levels eq:L2D.k/1989"),5]=-1The parameters PK1,...,PK5 also appear in the moment conditions related to y(i,t-1):
. mat gk[colnumb(e(Z), "Levels eq:L2D.y/1986"),1]=-bsys2[1,5]. mat gk[colnumb(e(Z), "Levels eq:L2D.y/1987"),1]=-bsys2[1,1]*bsys2[1,5]. mat gk[colnumb(e(Z), "Levels eq:L2D.y/1988"),1]=-bsys2[1,1]^2*bsys2[1,5]. mat gk[colnumb(e(Z), "Levels eq:L2D.y/1989"),1]=-bsys2[1,1]^3*bsys2[1,5] . mat gk[colnumb(e(Z), "Levels eq:L2D.y/1986"),2]=-bsys2[1,4]. mat gk[colnumb(e(Z), "Levels eq:L2D.y/1987"),2]=-(bsys2[1,1]*bsys2[1,4]+bsys2[1,5]). mat gk[colnumb(e(Z), "Levels eq:L2D.y/1988"),2]=-bsys2[1,1]*(bsys2[1,1]*bsys2[1,4]+bsys2[1,5]). mat gk[colnumb(e(Z), "Levels eq:L2D.y/1989"),2]=-bsys2[1,1]^2*(bsys2[1,1]*bsys2[1,4]+bsys2[1,5]) . mat gk[colnumb(e(Z), "Levels eq:L2D.y/1987"),3]=-bsys2[1,4]. mat gk[colnumb(e(Z), "Levels eq:L2D.y/1988"),3]=-(bsys2[1,1]*bsys2[1,4]+bsys2[1,5]). mat gk[colnumb(e(Z), "Levels eq:L2D.y/1989"),3]=-bsys2[1,1]*(bsys2[1,1]*bsys2[1,4]+bsys2[1,5]) . mat gk[colnumb(e(Z), "Levels eq:L2D.y/1988"),4]=-bsys2[1,4] . mat gk[colnumb(e(Z), "Levels eq:L2D.y/1989"),4]=-(bsys2[1,1]*bsys2[1,4]+bsys2[1,5]) . mat gk[colnumb(e(Z), "Levels eq:L2D.y/1989"),5]=-bsys2[1,4]. mat G=(G,gpsi,gn,gk)The test can be computed using the following formula:
. mat testcm = e(Ze)'*e(A2)*G*invsym(G'*e(A2)*G)*G'*e(A2)*e(Ze)We obtain a value of 33.32 with a p-value of 0.00046, so that we are led to reject the null hypothesis of validity of the mean stationarity assumption.