*! xtdolshm version 2.0
*! Performs panel data cointegration regression
*! Diallo Ibrahima Amadou
*! All comments are welcome, 16May2015






/* MAIN PROGRAM */

capture program drop xtdolshm
program xtdolshm, eclass byable(recall) sortpreserve
            version 10
            local options "Level(cilevel)"
            if replay() {
                         if "`e(cmd)'" != "xtdolshm" {
                                                      error 301
                         }
                         syntax [, `options']
            }
            else {
                  syntax varlist(ts) [if] [in] [, NLAgs(integer 2) NLEads(integer 1) `options' ]
                  marksample touse
                  qui count if `touse'
                  if r(N) == 0 {
                                di as err "No observations."
                                exit 2000
                  }
                  tempvar count
                  qui {
                       tsset
                       local ivar "`r(panelvar)'"
                       if "`ivar'" == "" {
                                          di as err "You must tsset the data: specify the panel and time variables."
                                          exit 459
                       }
                       bysort `touse' `ivar': gen long `count' = _N
                       sum `count' if `touse', meanonly
                       tempvar  mc
                       gen `mc' = .
                       replace `mc' = 1 if `count' == r(max) & `touse'
                       markout `touse' `mc'
                  }
                  gettoken first rest : varlist
                  qui tsset
                  _rmcoll `rest' if `touse'
                  local rest `r(varlist)'
                  qui tsset
                  local nrlags "`nlags'"
                  local nrleads "`nleads'"
                  ereturn clear
                  mata: _mz_dolsrd("`first'", "`rest'", "`r(panelvar)'","`touse'",`nrlags',`nrleads')
            }
            local sch = cond(!replay() , 0, 1)
            display
            DisplayDOLS, level(`level') sch(`sch')
            ereturn local cmdline "xtdolshm `0'"
end






/* DISPLAYS THE RESULTS */

program define DisplayDOLS, eclass sortpreserve
        syntax [, Level(cilevel) SCHwitch(integer 0) ]

        tempname bhat vce nng rr2 rr2a nlag nlead W  pvW omega12
        tempvar Ti touse

        local le "`level'"
        if `schwitch' == 1 {
                            matrix define `bhat' = e(b)
                            matrix define `vce' = e(V)
                            gen `touse' = e(sample)
        }
        else {
              matrix define `bhat' = e(bb)
              matrix define `vce' = e(VV)
              local touse = e(touse)
        }
        local nn = e(N)
        scalar `nng' = e(N_g)
        scalar `rr2' = e(r2)
        scalar `rr2a' = e(r2_a)
        local yy "`e(depvar)'"
        local ivar "`e(ivar)'"
        local k1 = e(k1)
        local T = e(T)
        scalar `nlag' =  e(nlags)
        scalar `nlead' = e(nleads)
        scalar `W' = e(chi2)
        scalar `pvW' = e(chi2_p)
        matrix define `omega12' = e(omega_1_2)
        ereturn clear
        sort `ivar'
        qui {
             preserve
             keep if `touse'
             by `ivar': gen long `Ti' = _N if _n ==_N
             summ `Ti'
             restore
        }
        local g_min "`r(min)'"
        local g_avg "`r(mean)'"
        local g_max "`r(max)'"
        ereturn clear
        ereturn post `bhat' `vce', esample(`touse')
        di _n in gr "DOLS Hom. Panel data Coint. Estimation results"                   _col(49) in gr "Number of obs" _col(68) "=" _col(70) in ye %9.0f `nn'
        di in gr "Group variable: " in ye abbrev("`ivar'",12)                         in gr _col(49) "Number of groups" _col(68) "=" _col(70) in ye %9.0g `nng'
        di in gr "Wald chi2(" as res `k1' as txt ")" _col(15) "= " as res %9.2f `W'   in gr _col(49) "Obs per group: min = " _col(70) in ye %9.0g `g_min'
        di in gr "Prob > chi2" _col(15) "=" as res %10.3f `pvW'                       in gr _col(49) "               avg = " _col(70) in ye %9.0g `g_avg'
        di in gr _col(49) "               max = " _col(70) in ye %9.0g `g_max'
        di in gr _col(49) "R-squared          =    " in ye %6.4f `rr2'
        di in gr _col(49) "Adj R-squared      =   " in ye %6.4f `rr2a'
        ereturn display, level(`le')
        ereturn matrix omega_1_2 `omega12'
        ereturn local depvar "`yy'"
        ereturn local ivar "`ivar'"
        ereturn scalar N = `nn'
        ereturn scalar k1 = `k1'
        ereturn scalar T = `T'
        ereturn scalar N_g = `nng'
        ereturn scalar r2 = `rr2'
        ereturn scalar r2_a = `rr2a'
        ereturn scalar g_avg = `g_avg'
        ereturn scalar g_min = `g_min'
        ereturn scalar g_max = `g_max'
        ereturn scalar nlags =  `nlag'
        ereturn scalar nleads = `nlead'
        ereturn scalar chi2 = `W'
        ereturn scalar chi2_p = `pvW'
        ereturn local predict xtdolshm_p
        ereturn local cmd "xtdolshm"
end






/*************************** PERFORMS THE DOLS ESTIMATIONS *****************************/
// This program is a MATA adaptation of the GAUSS original program of
// Chiang, M-H. and C. Kao(2002) " Nonstationary Panel Time Series Using NPT 1.3 -
// A User Guide," Center for Policy Research, Syracuse University.
// To keep the presentation simple and comparable to the original program, the notations
// are left as in the original program, which is downloadable at:
// http://faculty.maxwell.syr.edu/cdkao/working/npt.html
/***************************************************************************************/

mata:

void _mz_dolsrd(string scalar first, string scalar rest, string scalar panelvar,string scalar touse, real scalar nblags, real scalar nbleads)
{
              real matrix x0, x, y01, x01, info, temp_x1, x1, dif_x1, Dif_x1_cum, temp_x1_dols, x1_dols, temp_dif_x1, vb1, es1, ts1, /*
                      */ temp_D_resi_1, Dif_resi_1, /*
                      */ uk, ke, ini1, temp_sigma, ini2, amat_1, sigma_1, lr_var_1, epsilon_1, delta_1, /*
                      */ del_eps_1, DOLS_cov, R, r, num, tse, omega1_2, DOLS_cov_bon
              real colvector y0, pvar, y, b1, u1, tb1, DOLS_beta, u1_dols, y_dols
              real rowvector s, del_mu_eps_1, mu_eps_1
              string rowvector lvx
              string scalar lvy, lvp
              real scalar N, T, k1, lags, leads, Rsquare_1, Adjust_Rsquare_1, tss1, ess1, mu_1, mu_dot_eps_1, panel_kernel_lag_hom, /*
                         */  dof, ntot, W, pvW, scal1, DOLS_t_calc
              string matrix sg1



              lvy = tokens(first); lvx = tokens(rest); lvp = tokens(panelvar);
              st_view(y0,., lvy, touse);  st_view(x0,., st_tsrevar(lvx), touse); st_view(pvar,., lvp, touse)
              y01 = (pvar,y0); x01 = (pvar,x0)
              info = panelsetup(y01, 1, 2, 0); s = panelstats(info)
              N = s[1]; T = s[4]; k1 = cols(x0)
              y = mz_panrowshape(y01, 1, "KC", 1)
              x = mz_panrowshape(x01, 1, "KC", 1)
              lags  = nblags
              leads = nbleads
              y_dols = vec(mz_trimrows(y,(lags + 1),leads))
              y = y - mz_reshape(mean(y), T, N); y = vec(y)
              temp_x1 = J((T*N),1, missingof(x))
              for (i=1; i<=k1; i++) {
                                     temp_x1 = (temp_x1,vec(x[|1,(1+(i-1)*N)\.,(i*N)|]- mz_reshape(mean(x[|1,(1+(i-1)*N)\.,(i*N)|]),T,N)))
              }
              x1 = temp_x1[.,(2.. cols(temp_x1))]


              Dif_x1_cum = J((N*(T-1-(leads+lags))),1,missingof(x))
              dif_x1 = x[(2::rows(x)),.] - mz_trimrows(mz_decalen(x,1),1,0)
              for (i=0; i<=(leads+lags); i++) {
                                               if (i != lags) {
                                                               for (j=1; j<=k1; j++) {
                                                                                      Dif_x1_cum = (Dif_x1_cum,vec(mz_trimrows(dif_x1[|1,(1+(j-1)*N)\.,(j*N)|],i,(leads + lags - i))))
                                                               }
                                               }
              }
              temp_x1_dols = J((N*(T-1-(leads+lags))),1,missingof(x))
              for (j=1; j<=k1; j++) {
                                     temp_x1_dols = (temp_x1_dols,vec(mz_trimrows(x[|1,(1+(j-1)*N)\.,(j*N)|],(lags+1),leads)))
              }
              x1_dols = (vec(J((T-(leads+lags+1)),N,1)),temp_x1_dols[.,(2.. cols(temp_x1_dols))],Dif_x1_cum[.,(2.. cols(Dif_x1_cum))])

              temp_dif_x1 = J((N*(T-1)),1,missingof(x))
              for (i=1; i<=k1; i++) {
                                     temp_dif_x1 = (temp_dif_x1,vec(dif_x1[|1,(1+(i-1)*N)\.,(i*N)|]))
              }
              dif_x1 = temp_dif_x1[.,(2.. cols(temp_dif_x1))]
              b1 = invsym(cross(x1, x1))*cross(x1, y)
              u1 = y - x1*b1
              vb1 = invsym(cross(x1, x1))*(cross(u1,u1)/(N*T-k1))
              tb1 = b1:/diagonal(sqrt(vb1))
              es1 = mz_reshape(u1,N,T)'; ess1 = colsum((T-1)*(mz_stdc(es1):*mz_stdc(es1)))'
              ts1 = mz_reshape(y,N,T)';  tss1 = colsum((T-1)*(mz_stdc(ts1):*mz_stdc(ts1)))'
              Rsquare_1 = 1-ess1/tss1
              Adjust_Rsquare_1 = 1-((cross(u1,u1)/(N*T-k1))/(tss1/(N*T)))
              temp_D_resi_1 = mz_reshape(u1,N,T)
              Dif_resi_1 = temp_D_resi_1[.,(2.. T)]'
              ini1 = ini2 = amat_1 = sigma_1 = lr_var_1 /*
              */ = J((cols(dif_x1) + 1),(cols(dif_x1) + 1),0)
              panel_kernel_lag_hom = 5

              for (i=1; i<=N; i++) {
                                    uk = (Dif_resi_1[.,i],dif_x1[|(1+(i-1)*(T-1)),1\(i*(T-1)),.|])
                                    ke = mz_fejer(uk, panel_kernel_lag_hom)
                                    ini1 = ini1 + ke
                                    temp_sigma = cross(uk,uk)/rows(uk)
                                    ini2 = ini2 + temp_sigma
              }

              amat_1 = ini1/N
              sigma_1 = ini2/N
              lr_var_1 = sigma_1 + amat_1 + amat_1'
              epsilon_1 = lr_var_1[(2::rows(lr_var_1)),(2.. cols(lr_var_1))]
              delta_1 = sigma_1 + amat_1
              del_mu_eps_1 = delta_1[1,(2.. cols(delta_1))]
              del_eps_1 = delta_1[(2::rows(delta_1)),(2.. cols(delta_1))]
              mu_1 = lr_var_1[1,1]
              mu_eps_1 = lr_var_1[1,(2.. cols(lr_var_1))]
              mu_dot_eps_1 = mu_1 - mu_eps_1*qrinv(epsilon_1)*mu_eps_1'
              DOLS_beta = invsym(cross(x1_dols, x1_dols))*cross(x1_dols, y_dols)
              DOLS_cov = 6*qrinv(epsilon_1)*mu_dot_eps_1
              DOLS_cov_bon = DOLS_cov/((sqrt(N)*(T - (leads + lags + 1)))^2)
              u1 = y_dols - x1_dols*DOLS_beta
              u1_dols = y - x1*DOLS_beta[|2\(k1+1)|]
              es1 = mz_reshape((x1_dols*DOLS_beta),N,(T-(leads + lags + 1)))'
              ess1 = colsum((T - (leads + lags + 2))*(mz_stdc(es1) :* mz_stdc(es1)))'
              ts1 = mz_reshape(y_dols,N,(T - (leads + lags + 1)))'
              tss1 = colsum((T - (leads + lags + 2))*(mz_stdc(ts1) :* mz_stdc(ts1)))'
              Rsquare_1 = ess1/tss1
              Adjust_Rsquare_1 = 1 - /*
                   */ ((cross(u1,u1)/(N*(T - (leads + lags + 1)) - N - k1 - (lags + leads )*k1))/ /*
                   */ (tss1/(N*(T - (leads + lags + 1)) - N)))
              omega1_2 = mu_dot_eps_1
              ntot = N*(T-(leads + lags + 1))
              DOLS_betaf = DOLS_beta[|2\(k1+1)|]'
              R = I(k1)
              r = J(k1,1,0)
              W = (1/6)*N*((T - (leads + lags + 1))^2)*((R*DOLS_betaf'-r)')*qrinv(R*qrinv(epsilon_1)*mu_dot_eps_1*R')*(R*DOLS_betaf'-r)
              pvW = chi2tail(k1, W)
              st_matrix("e(bb)", DOLS_betaf)
              st_matrix("e(VV)", DOLS_cov_bon)
              st_matrix("e(omega_1_2)",omega1_2)
              st_global("e(touse)", touse)
              sg1 = ("")
              sg2 = (mz_reshape(sg1,k1,1), lvx')
              st_matrixcolstripe("e(bb)", sg2)
              st_matrixcolstripe("e(VV)", sg2)
              st_matrixrowstripe("e(VV)", sg2)
              st_numscalar("e(N)", ntot)
              st_numscalar("e(N_g)", N)
              st_numscalar("e(r2)", Rsquare_1)
              st_numscalar("e(r2_a)", Adjust_Rsquare_1)
              st_global("e(depvar)", lvy)
              st_global("e(ivar)", lvp)
              st_numscalar("e(T)", T)
              st_numscalar("e(k1)", k1)
              st_numscalar("e(nlags)", lags)
              st_numscalar("e(nleads)", leads)
              st_numscalar("e(chi2)", W)
              st_numscalar("e(chi2_p)", pvW)
}


real matrix mz_fejer(real matrix x, real scalar k)
{
     real matrix t1, t2, a
     real scalar i, m, f
     a = J(cols(x),cols(x),0)
     for(i=1; i<=k; i++) {
                          f = i/(k+1)
                          m = 1 - f
                          t1 = mz_trimrows(x,i,0)
                          t2 = mz_trimrows(mz_decalen(x,i),i,0)
                          a = a + m*cross(t1,t2)
     }
     return(a/rows(x))
}



end