Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • - Acreg - spatial use: Missing standard errors for certain distance cutoffs

    Dear all,

    I encounter the problem with the user written command acreg or ols_spatial_HAC for spatial standard error adjustments that standard errors are not displayed for certain ranges of distance cutoffs. I use STATA 15 on Win (64-bit).

    I have 500+ observations but this reduced example should show the problem:

    Code:
    *My example data
    clear all
    input float(lon lat) byte x double y float constant
     8.978752 49.15352 1             .3125 1
      9.16124 47.79682 0  .230769230769231 1
     8.627637 48.93355 1  .277777777777778 1
     9.386187 47.71273 0  .333333333333333 1
     9.283757 49.29697 1 .0952380952380952 1
     7.943795 48.47777 1                .3 1
     8.503139 49.29774 1  .285714285714286 1
     8.572759 48.94195 1  .272727272727273 1
     7.883648 48.07095 1  .227272727272727 1
     8.195633 47.67924 0  .166666666666667 1
     7.950767 47.96292 0  .333333333333333 1
     8.760374 47.78041 1                .2 1
     8.185254 48.86308 1  .234042553191489 1
    8.3561125 48.71901 1  .166666666666667 1
     8.644542 49.05665 1  .214285714285714 1
     7.594427 47.63269 1                .5 1
     8.70044 49.34214 1  .307692307692308 1
    end
    *now run my regression for various distance cutoff values
    forvalues i=25(25) 200{
    acreg y x , spatial latitude(lat) longitude(lon) distcutoff(`i')
    }
    *I cannot explain now why there is no se shown for 175 but one for 200 exists
    
    *by adding the bartlett option this issue disappers
    forvalues i=25(25) 200{
    acreg y x , spatial latitude(lat) longitude(lon) distcutoff(`i') bartlett
    }
    This also happens when using ols_spatial_HAC command instead of acreg.

    I am puzzled by the fact that I get a standard error for the 200km distance cutoff, but none for the 175km one.


    Thank you in advance!

    Sincerely,
    Steven
    Last edited by Steven Gregory Hill; 22 Oct 2019, 16:57.

  • #2
    Hi Steve,
    You should know that acreg command is hardly used by the spatial modeler, it is new and developed by one PhD student at swiss Uni. I have seen PhD student from Monash University using this code. I don't think you will get any useful reply here, I suggest you to contact the author.

    Comment


    • #3
      On this occasion I doubt that example data will clarify. My guess is that this is just a side-effect of the spatial distribution of the data. You don't have enough information in that distance range. Just a guess, but it's all you have so far,

      Comment


      • #4
        Hi Harry,

        thank you. I also tried this with the Code by Hsiang (PNAS 2010) with the same result.
        Maybe more people are using this user-written command?

        Code:
        clear
        input float(lon lat) byte x double y float(constant year panel touse) byte(_est_OLS _est_spatial _est_spatHAC)
         8.67651 48.94452 1  .285714285714286 1 1  1 1 1 1 1
        8.989747 49.20008 1 .0666666666666667 1 1  2 1 1 1 1
         8.90716 49.13757 1  .222222222222222 1 1  3 1 1 1 1
        7.956594 47.56825 0  .208333333333333 1 1  4 1 1 1 1
        7.776893 48.06351 1  .111111111111111 1 1  5 1 1 1 1
        8.223049 48.35345 0  .416666666666667 1 1  6 1 1 1 1
        7.700581 47.58864 1  .166666666666667 1 1  7 1 1 1 1
        9.035539 47.81129 0             .1875 1 1  8 1 1 1 1
        8.401457 47.93015 0                .1 1 1  9 1 1 1 1
        9.196566 49.46891 0  .130434782608696 1 1 10 1 1 1 1
        8.694172 49.40551 1             .3125 1 1 11 1 1 1 1
        8.536718 49.14109 1  .333333333333333 1 1 12 1 1 1 1
        9.324537 49.43767 1              .125 1 1 13 1 1 1 1
        8.078202 47.89349 0  .166666666666667 1 1 14 1 1 1 1
        8.754618 48.96887 1  .333333333333333 1 1 15 1 1 1 1
         7.77809 47.64332 1  .142857142857143 1 1 16 1 1 1 1
        8.506112 49.39533 1  .227272727272727 1 1 17 1 1 1 1
        8.097216 47.78621 0  .166666666666667 1 1 18 1 1 1 1
        9.427267 49.58295 0  .133333333333333 1 1 19 1 1 1 1
        7.769179 47.58644 1              .375 1 1 20 1 1 1 1
         8.00207 48.71616 1              .125 1 1 21 1 1 1 1
        8.475474 49.09151 1  .307692307692308 1 1 22 1 1 1 1
        8.535838 49.05559 1  .277777777777778 1 1 23 1 1 1 1
        8.601546 49.02109 1  .333333333333333 1 1 24 1 1 1 1
        7.979859 48.55326 1  .277777777777778 1 1 25 1 1 1 1
         9.07665  49.3006 1  .230769230769231 1 1 26 1 1 1 1
        8.522661 49.29453 1  .388888888888889 1 1 27 1 1 1 1
        7.807001 48.12419 1  .206896551724138 1 1 28 1 1 1 1
        7.849948 48.57199 1  .192307692307692 1 1 29 1 1 1 1
        8.418171 47.63171 0  .181818181818182 1 1 30 1 1 1 1
        7.886527 48.03259 1  .409090909090909 1 1 31 1 1 1 1
        7.918152 48.11533 0  .333333333333333 1 1 32 1 1 1 1
        8.503129 48.90504 1               .28 1 1 33 1 1 1 1
        7.683842 47.62807 1            .40625 1 1 34 1 1 1 1
        9.069303 49.39884 1  .166666666666667 1 1 35 1 1 1 1
        8.479564 49.15499 1  .222222222222222 1 1 36 1 1 1 1
        7.733454  47.8698 1               .25 1 1 37 1 1 1 1
        8.296137 48.41405 0                .1 1 1 38 1 1 1 1
        7.812906 48.16362 1  .333333333333333 1 1 39 1 1 1 1
        7.987206  48.1443 0  .357142857142857 1 1 40 1 1 1 1
        9.635777 49.38785 0  .181818181818182 1 1 41 1 1 1 1
        8.799194 48.90755 1  .272727272727273 1 1 42 1 1 1 1
        7.818086 47.99253 1  .333333333333333 1 1 43 1 1 1 1
        8.884067 49.37471 1 .0666666666666667 1 1 44 1 1 1 1
        8.552239 49.19415 1  .285714285714286 1 1 45 1 1 1 1
        8.044605 48.63147 1  .269230769230769 1 1 46 1 1 1 1
        8.407392 49.12775 1  .166666666666667 1 1 47 1 1 1 1
         8.02769 47.96811 0               .25 1 1 48 1 1 1 1
        8.121563 47.62234 0  .107142857142857 1 1 49 1 1 1 1
        end
        
        forvalues i=25(25) 300{
        ols_spatial_HAC y x constant, lat(lat) lon(lon) timevar(year) panelvar(panel) distcutoff(`i')
        }
        What alternative code would you advise me to use instead?

        Comment


        • #5
          Please tell us where these community-contributed commands can be found. That's a longstanding request here. https://www.statalist.org/forums/help#stata 12.1 spells it out.

          Comment


          • #6
            Hello Nick,

            thank you for your reply. I am sorry, I was under the false impression that the hyperlink above would lead you to the code.
            I looked the code up on this website of the author, Solomon Hsiang, here: http://www.fight-entropy.com/2010/06...t-ols-for.html


            The ado file he has produced on which the command ols_spaital_HAC is build is the following:

            Code:
             
            *! S. HSIANG 6/2010: PROGRAM TO ESTIMATE SPATIAL HAC ERRORS FOR OLS REGRESSION MODEL [V3 UPDATE 6/2018]
            
            /*-----------------------------------------------------------------------------
            
             v1 S. HSIANG 6/10 [[email protected]]
            
             v2 UPDATE 6/13 [[email protected]]:
                INTRODUCED 'DROPVAR' OPTION BASED ON CODE PROVIDED BY KYLE MENG
             
             V3 UPDATE 6/18 [[email protected]]:
                CORRECTED ERROR (LINE 428) FOUND BY MATHIAS THOENIG THAT INCORRECTLY
                COMPUTED WEIGHTS FOR INTER-TEMPORAL AUTOCORRELATION ESTIMATES WITHIN PANEL
                UNITS. SIGN OF BIAS IN V2 IS INDETERMINATE, DEPENDS ON LAG LENGTH AND DATA
                STRUCTURE.
             
             ------------------------------------------------------------------------------
            
             This may contain errors. Please notify me of any errors you find.
             
             ------------------------------------------------------------------------------
            
             Syntax:
             
             ols_spatial_HAC Yvar Xvarlist, lat(latvar) lon(lonvar) Timevar(tvar) Panelvar(pvar) [DISTcutoff(#) LAGcutoff(#) bartlett DISPlay star dropvar]
            
             Function calculates non-parametric (GMM) spatial and autocorrelation
             structure using a panel data set.  Spatial correlation is estimated for all
             observations within a given period.  Autocorrelation is estimated for a
             given individual over multiple periods up to some lag length. Var-Covar
             matrix is robust to heteroskedasticity.
             
             A variable equal to 1 is required to estimate a constant term.
             
             Example commands:
             
             ols_spatial_HAC dep indep1 indep2 const, lat(C1) lon(C2) t(year) p(id) dist(300) lag(3) bartlett disp
            
             ols_spatial_HAC dep indep*, lat(C1) lon(C2) timevar(year) panelvar(id) dist(100) lag(2) star dropvar
            
             ------------------------------------------------------------------------------
             
             Requred arguments:
             
             Yvar: dependent variable  
             Xvarlist: independnet variables (INCLUDE constant as column)
             latvar: variable containing latitude in DEGREES of each obs
             lonvar: same, but longitude
             tvar: varible containing time variable
             pvar: variable containing panel variable (must be numeric, see "encode")
             
             ------------------------------------------------------------------------------
             
             Optional arguments:
             
             distcutoff(#): {abbrev dist(#)} describes the distance cutoff in KILOMETERS for the spatial kernal (the distance at which spatial correlation is assumed to vanish). Default is 1 KM.
             
             lagcutoff(#): {abbrev lag(#)} describes the maximum number of temporal periods for the linear Bartlett window that weights serial correlation across time periods (the distance at which serial correlation is assumed to vanish). Default is 0 PERIODS (no serial correlation). {Note, Greene recommends at least T^0.25}  
             
             ------------------------------------------------------------------------------
             
             Options:
             
             bartlett: use a linear bartlett window for spatial correlations, instead of a uniform kernal
             
             display: {abbrev disp} display a table with estimated coeff and SE & t-stat using OLS, adjusting for spatial correlation and adjusting for both spatial and serial correlation. Can be used with star option. Ex:
             
             -----------------------------------------------
                 Variable |   OLS      spatial    spatHAC   
             -------------+---------------------------------
                   indep1 |    0.568      0.568      0.568  
                          |    0.198      0.206      0.240  
                          |    2.876      2.761      2.369  
                    const |    6.415      6.415      6.415  
                          |    0.790      1.176      1.340  
                          |    8.119      5.454      4.786  
             -----------------------------------------------
                                              legend: b/se/t
             
            
             star: same as display, but uses stars to denote significance and does not show SE & t-stat. Can be used with display option. Ex:
             
             -----------------------------------------------------
                 Variable |    OLS        spatial      spatHAC    
             -------------+---------------------------------------
                   indep1 |   0.568***     0.568***     0.568**   
                    const |   6.415***     6.415***     6.415***  
             -----------------------------------------------------
                               legend: * p<.1; ** p<.05; *** p<.01
                               
                               
             dropvar: Drops variables that Stata would drop due to collinearity. This requires that an additiona regression is run, so it slows the code down. For large datasets, if this function is called many times, it may be faster to ensure that colinear variables are dropped in advance rather than using the option dropvar. If Stata returns "estimates post: matrix has missing values", than including the option dropvar may solve the problem. (This option written by Kyle Meng).
             
             ------------------------------------------------------------------------------
             
             Implementation:
             
             The default kernal used to weight spatial correlations is a uniform kernal that
             discontinously falls from 1 to zero at length locCutoff in all directions (it is isotropic). This is the kernal recommented by Conley (2008). If the option "bartlett" is selected, a conical kernal that decays linearly with distance in all directions is used instead.
             
             Serial correlation bewteen observations of the same individual over multiple periods seperated by lag L are weighted by
            
                   w(L) = 1 - L/(lagCutoff+1)
                   
             ------------------------------------------------------------------------------
            
             Notes:
            
             Location arguments should specify lat-lon units in DEGREES, however
             distcutoff should be specified in KILOMETERS.
            
             distcutoff must exceed zero. CAREFUL: do not supply
             coordinate locations in modulo(360) if observations straddle the
             zero-meridian or in modulo(180) if they straddle the date-line.
            
             Distances are computed by approximating the planet's surface as a plane
             around each observation.  This allows for large changes in LAT to be
             present in the dataset (it corrects for changes in the length of
             LON-degrees associated with changes in LAT). However, it does not account
             for the local curvature of the surface around a point, so distances will
             be slightly lower than true geodesics. This should not be a concern so
             long as locCutoff is < O(~2000km), probably.
            
             Each time-series for an individual observation in the panel is treated
             with Heteroskedastic and Autocorrelation Standard Errors. If lagcutoff =
             0, than this estimate is equivelent to White standard errors (with spatial correlations
             accounted for). If lagcutoff = infinity, than this treatment is
             equivelent to the "cluster" command in Stata at the panel variable level.
            
             This script stores estimation results in standard Stata formats, so most "ereturn" commands should work properly.  It is also compatible with "outreg2," although I have not tested other programs.
            
             The R^2 statistics output by this function will differ from analogous R^2 stats
             computed using "reg" since this function omits the constant.
             ------------------------------------------------------------------------------
            
             References:
            
                  TG Conley "GMM Estimation with Cross Sectional Dependence"
                  Journal of Econometrics, Vol. 92 Issue 1(September 1999) 1-45
                  http://www.elsevier.com/homepage/sae/econworld/econbase/econom/frame.htm
                  
                  and
            
                  Conley "Spatial Econometrics" New Palgrave Dictionary of Economics,
                  2nd Edition, 2008
            
                  and
            
                  Greene, Econometric Analysis, p. 546
            
                  and
            
                  Modified from scripts written by Ruben Lebowski and Wolfram Schlenker and Jean-Pierre Dube and Solomon Hsiang
                  Debugging help provided by Mathias Thoenig.
             
             -----------------------------------------------------------------------------*/
            
            program ols_spatial_HAC, eclass byable(recall)
            version 11
            syntax varlist(ts fv min=2) [if] [in], ///
                            lat(varname numeric) lon(varname numeric) ///
                            Timevar(varname numeric) Panelvar(varname numeric) [LAGcutoff(integer 0) DISTcutoff(real 1) ///
                            DISPlay star bartlett dropvar]
                            
            /*--------PARSING COMMANDS AND SETUP-------*/
            
            capture drop touse
            marksample touse                // indicator for inclusion in the sample
            gen touse = `touse'
            
            //parsing variables
            loc Y = word("`varlist'",1)        
            
            loc listing "`varlist'"
            
            loc X ""
            scalar k = 0
            
            //make sure that Y is not included in the other_var list
            foreach i of loc listing {
                if "`i'" ~= "`Y'"{
                    loc X "`X' `i'"
                    scalar k = k + 1 // # indep variables
                    
                }
            }
            
            
            //Kyle Meng's code to drop omitted variables that Stata would drop due to collinearity
            
            if "`dropvar'" == "dropvar"{
                
                quietly reg `Y' `X' if `touse', nocons
                
                mat omittedMat=e(b)
                local newVarList=""
                local i=1
                scalar k = 0 //replace the old k if this option is selected
                
                foreach var of varlist `X'{
                    if omittedMat[1,`i']!=0{
                        loc newVarList "`newVarList' `var'"
                        scalar k = k + 1
                    }
                    local i=`i'+1
                }
                
                loc X "`newVarList'"
            }
            
            //generating a function of the included obs
            quietly count if `touse'        
            scalar n = r(N)                    // # obs
            scalar n_obs = r(N)
            
            /*--------FIRST DO OLS, STORE RESULTS-------*/
            
            
            quietly: reg `Y' `X' if `touse', nocons
            estimates store OLS
            
            //est tab OLS, stats(N r2)
            
            /*--------SECOND, IMPORT ALL VALUES INTO MATA-------*/
            
            mata{
            
            Y_var = st_local("Y") //importing variable assignments to mata
            X_var = st_local("X")
            lat_var = st_local("lat")
            lon_var = st_local("lon")
            time_var = st_local("timevar")
            panel_var = st_local("panelvar")
            
            //NOTE: values are all imported as "views" instead of being copied and pasted as Mata data because it is faster, however none of the matrices are changed in any way, so it should not permanently affect the data.
            
            st_view(Y=.,.,tokens(Y_var),"touse") //importing variables vectors to mata
            st_view(X=.,.,tokens(X_var),"touse")
            st_view(lat=.,.,tokens(lat_var),"touse")
            st_view(lon=.,.,tokens(lon_var),"touse")
            st_view(time=.,.,tokens(time_var),"touse")
            st_view(panel=.,.,tokens(panel_var),"touse")
            
            k = st_numscalar("k")                //importing other parameters
            n = st_numscalar("n")
            b = st_matrix("e(b)")                // (estimated coefficients, row vector)
            lag_var = st_local("lagcutoff")
            lag_cutoff = strtoreal(lag_var)
            dist_var = st_local("distcutoff")
            dist_cutoff = strtoreal(dist_var)
            
            XeeX = J(k, k, 0)                 //set variance-covariance matrix equal to zeros
            
            
            /*--------THIRD, CORRECT VCE FOR SPATIAL CORR-------*/
            
            timeUnique = uniqrows(time)
            Ntime = rows(timeUnique)         // # of obs. periods
            
            for (ti = 1; ti <= Ntime; ti++){
                
                
            
                // 1 if in year ti, 0 otherwise:
            
                rows_ti = time:==timeUnique[ti,1]     
            
                //get subsets of variables for time ti (without changing original matrix)
                
                Y1 = select(Y, rows_ti)
                X1 = select(X, rows_ti)
                lat1 = select(lat, rows_ti)
                lon1 = select(lon, rows_ti)
                e1 = Y1 - X1*b'
                
                n1 = length(Y1)             // # obs for period ti
                
                //loop over all observations in period ti
            
                for (i = 1; i <=n1; i++){
                    
            
                    //----------------------------------------------------------------
                    // step a: get non-parametric weight
                
                    //This is a Euclidean distance scale IN KILOMETERS specific to i
                    
                    lon_scale = cos(lat1[i,1]*pi()/180)*111
                    lat_scale = 111
                    
            
                    // Distance scales lat and lon degrees differently depending on
                    // latitude.  The distance here assumes a distortion of Euclidean
                    // space around the location of 'i' that is approximately correct for
                    // displacements around the location of 'i'
                    //
                    //    Note:     1 deg lat = 111 km
                    //             1 deg lon = 111 km * cos(lat)
                    
                    distance_i = ((lat_scale*(lat1[i,1]:-lat1)):^2 + ///     
                                  (lon_scale*(lon1[i,1]:-lon1)):^2):^0.5
            
            
                    
                    // this sets all observations beyon dist_cutoff to zero, and weights all nearby observations equally [this kernal is isotropic]
                    
                    window_i = distance_i :<= dist_cutoff
            
                    //----------------------------------------------------------------
                    // adjustment for the weights if a "bartlett" kernal is selected as an option
             
                    if ("`bartlett'"=="bartlett"){
                    
                        // this weights observations as a linear function of distance
                        // that is zero at the cutoff distance
                        
                        weight_i = 1:- distance_i:/dist_cutoff
            
                        window_i = window_i:*weight_i
                    }
            
             
                    //----------------------------------------------------------------
                    // step b: construct X'e'eX for the given observation
             
                     XeeXh = ((X1[i,.]'*J(1,n1,1)*e1[i,1]):*(J(k,1,1)*e1':*window_i'))*X1
            
                    //add each new k x k matrix onto the existing matrix (will be symmetric)
                
                    XeeX = XeeX + XeeXh     
                
                } //i
            } // ti
            
            
            
            // -----------------------------------------------------------------
            // generate the VCE for only cross-sectional spatial correlation,
            // return it for comparison
            
            invXX = luinv(X'*X) * n
            
            XeeX_spatial = XeeX / n
            
            V = invXX * XeeX_spatial * invXX / n
            
            // Ensures that the matrix is symmetric
            // in theory, it should be already, but it may not be due to rounding errors for large datasets
            V = (V+V')/2
            
            st_matrix("V_spatial", V)
            
            } // mata
            
            
            //------------------------------------------------------------------
            // storing old statistics about the estimate so postestimation can be used
            
            matrix beta = e(b)
            scalar r2_old = e(r2)
            scalar df_m_old = e(df_m)
            scalar df_r_old = e(df_r)
            scalar rmse_old = e(rmse)
            scalar mss_old = e(mss)
            scalar rss_old = e(rss)
            scalar r2_a_old = e(r2_a)
            
            // the row and column names of the new VCE must match the vector b
            
            matrix colnames V_spatial = `X'
            matrix rownames V_spatial = `X'
             
            // this sets the new estimates as the most recent model
            
            ereturn post beta V_spatial, esample(`touse')
            
            // then filling back in all the parameters for postestimation
            
            ereturn local cmd = "ols_spatial"
            
            ereturn scalar N = n_obs
            
            ereturn scalar r2 = r2_old
            ereturn scalar df_m = df_m_old
            ereturn scalar df_r = df_r_old
            ereturn scalar rmse = rmse_old
            ereturn scalar mss = mss_old
            ereturn scalar rss = rss_old
            ereturn scalar r2_a = r2_a_old
            
            ereturn local title = "Linear regression"
            ereturn local depvar = "`Y'"
            ereturn local predict = "regres_p"
            ereturn local model = "ols"
            ereturn local estat_cmd = "regress_estat"
            
            //storing these estimates for comparison to OLS and the HAC estimates
            
            estimates store spatial
            
            
            
            /*--------FOURTH, CORRECT VCE FOR SERIAL CORR-------*/
            
            mata{
            
            panelUnique = uniqrows(panel)
            Npanel = rows(panelUnique)         // # of panels
            
            for (pi = 1; pi <= Npanel; pi++){
                
                // 1 if in panel pi, 0 otherwise:
            
                rows_pi = panel:==panelUnique[pi,1]     
            
                //get subsets of variables for panel pi (without changing original matrix)
                
                Y1 = select(Y, rows_pi)
                X1 = select(X, rows_pi)
                time1 = select(time, rows_pi)
                e1 = Y1 - X1*b'
            
                n1 = length(Y1)             // # obs for panel pi
                
                //loop over all observations in panel pi
            
                for (t = 1; t <=n1; t++){
            
                       // ----------------------------------------------------------------
                    // step a: get non-parametric weight
                    
                    // this is the weight for Newey-West with a Bartlett kernal
                    
                    //weight = (1:-abs(time1[t,1] :- time1))/(lag_cutoff+1) // correction: need to removing parentheses to compute inter-temporal  (6/10/18)
                    weight = 1:-abs(time1[t,1] :- time1)/(lag_cutoff+1)
            
                    
                    // obs var far enough apart in time are prescribed to have no estimated
                    // correlation (Greene recomments lag_cutoff >= T^0.25 {pg 546})
                    
                    window_t = (abs(time1[t,1]:- time1) :<= lag_cutoff) :* weight
                    
                    //this is required so diagonal terms in var-covar matrix are not
                    //double counted (since they were counted once above for the spatial
                    //correlation estimates:
                    
                    window_t = window_t :* (time1[t,1] :!= time1)                   
                        
                      // ----------------------------------------------------------------
                    // step b: construct X'e'eX for given observation
                     
                       XeeXh = ((X1[t,.]'*J(1,n1,1)*e1[t,1]):*(J(k,1,1)*e1':*window_t'))*X1
            
                    //add each new k x k matrix onto the existing matrix (will be symmetric)
                            
                    XeeX = XeeX + XeeXh
            
                } // t
            } // pi
            
            
            
            
            // -----------------------------------------------------------------
            // generate the VCE for x-sectional spatial correlation and serial correlation
            
            XeeX_spatial_HAC = XeeX / n
            
            V = invXX * XeeX_spatial_HAC * invXX / n
            
            // Ensures that the matrix is symmetric
            // in theory, it should be already, but it may not be due to rounding errors for large datasets
            V = (V+V')/2
            
            st_matrix("V_spatial_HAC", V)
            
            } // mata
            
            //------------------------------------------------------------------
            //storing results
            
            matrix beta = e(b)
            
            // the row and column names of the new VCE must match the vector b
            
            matrix colnames V_spatial_HAC = `X'
            matrix rownames V_spatial_HAC = `X'
            
            // this sets the new estimates as the most recent model
            
            marksample touse                // indicator for inclusion in the sample
            
            ereturn post beta V_spatial_HAC, esample(`touse')
            
            // then filling back in all the parameters for postestimation
            
            ereturn local cmd = "ols_spatial_HAC"
            
            ereturn scalar N = n_obs
            ereturn scalar r2 = r2_old
            ereturn scalar df_m = df_m_old
            ereturn scalar df_r = df_r_old
            ereturn scalar rmse = rmse_old
            ereturn scalar mss = mss_old
            ereturn scalar rss = rss_old
            ereturn scalar r2_a = r2_a_old
            
            ereturn local title = "Linear regression"
            ereturn local depvar = "`Y'"
            ereturn local predict = "regres_p"
            ereturn local model = "ols"
            ereturn local estat_cmd = "regress_estat"
            
            //storing these estimates for comparison to OLS and the HAC estimates
            
            estimates store spatHAC
            
            //------------------------------------------------------------------
            //displaying results
            
            disp as txt " "
            disp as txt "OLS REGRESSION"
            disp as txt " "
            disp as txt "SE CORRECTED FOR CROSS-SECTIONAL SPATIAL DEPENDANCE"
            disp as txt "             AND PANEL-SPECIFIC SERIAL CORRELATION"
            disp as txt " "
            disp as txt "DEPENDANT VARIABLE: `Y'"
            disp as txt "INDEPENDANT VARIABLES: `X'"
            disp as txt " "
            disp as txt "SPATIAL CORRELATION KERNAL CUTOFF: `distcutoff' KM"
            
            if "`bartlett'" == "bartlett" {
                disp as txt "(NOTE: LINEAR BARTLETT WINDOW USED FOR SPATIAL KERNAL)"
            }
                
            disp as txt "SERIAL CORRELATION KERNAL CUTOFF: `lagcutoff' PERIODS"
            
            ereturn display // standard Stata regression table format
            
            // displaying different SE if option selected
            
            if "`display'" == "display"{
                disp as txt " "
                disp as txt "STANDARD ERRORS UNDER OLS, WITH SPATIAL CORRECTION AND WITH SPATIAL AND SERIAL CORRECTION:"
                estimates table OLS spatial spatHAC, b(%7.3f) se(%7.3f) t(%7.3f) stats(N r2)     
            }
            
            if "`star'" == "star"{
                disp as txt " "
                disp as txt "STANDARD ERRORS UNDER OLS, WITH SPATIAL CORRECTION AND WITH SPATIAL AND SERIAL CORRECTION:"
                estimates table OLS spatial spatHAC, b(%7.3f) star(0.10 0.05 0.01)
            }
            
            //------------------------------------------------------------------
            // cleaning up Mata environment
            
            capture mata mata drop V invXX  XeeX XeeXh XeeX_spatial_HAC window_t window_i weight t i ti pi X1 Y1 e1 time1 n1 lat lon lat1 lon1 lat_scale lon_scale rows_ti rows_pi timeUnique panelUnique Ntime Npanel X X_var XeeX_spatial Y Y_var b dist_cutoff dist_var distance_i k lag_cutoff lag_var lat_var lon_var n panel panel_var time time_var weight_i
            
            /*
            if "`bartlett'" == "bartlett" {
                capture mata mata drop weight_i            
            }
            */
            
            end

            Comment

            Working...
            X