*! version 0.9.0 hatemij

Title

hatemij -- Hatemi-J test for cointegration with two unknown breaks

Syntax

hatemij
depvar
indepvars
[if] [in]
[, options]

Main
model(string) deterministic specification: c = constant, ct = constant and trend, rs = regime shift; default model(rs)

choice(#) ADF lag selection: 1 = fixed (kmax), 2 = AIC, 3 = BIC, 4 = downward-t (1.96), 5 = downward-t (1.645), 6 = Breusch-Godfrey test based selection; default choice(2)

kmax(#) maximum lag for the ADF test; default kmax(12)

kernel(string) iid, bartlett, or qs; default kernel(qs)

bwl(#) fixed bandwidth for the long-run variance; default is the Andrews (1991) automatic bandwidth

trimming(#) trimming rate for the break search; default trimming(0.15)

bglags(#) Breusch-Godfrey test horizon

reg(string) reports cointegrating regression and add D1, D2 and interaction terms to the dataset

regnewey(string) as reg() but with Newey-West HAC standard errors (newey)

tspdlib compute the statistics with Gauss's TSPDLIB coint_hatemiJ routine instead of Hatemi-J's coint2b routine of Gauss

The data must be tsset. depvar and indepvars must be numeric.
indepvars may contain at most four variables: Hatemi-J (2008) tabulates
critical values only for up to four regressors (m <= 4), so hatemij stops
with an error if more than four are specified.

Description

The hatemij command implements all three models of the Hatemi-J
cointegration test. It also provides a model-selection algorithm that chooses a
specification free from autocorrelation [choice(6) option]. This approach was previously used in the
kapetanios command and can also be applied through hatemij. The
algorithm is highly efficient and typically produces results within a few
seconds.

hatemij implements the residual-based tests for cointegration with two
unknown regime shifts proposed by Hatemi-J (2008). The deterministic part is
chosen with model(). Writing x(t) for the regressors and D1(t), D2(t) for
the two level-shift dummies (each equal to 0 before its break date and 1
afterwards, with both break dates chosen from the data), the three
specifications are:

model(c) (constant): y(t) = a0 + a1 D1(t) + a2 D2(t) + b'x(t) + u(t)

model(ct) (constant and trend): y(t) = a0 + a1 D1(t) + a2 D2(t) + g t + b'x(t) + u(t)

model(rs) (regime shift): y(t) = a0 + a1 D1(t) + a2 D2(t) + b0'x(t) + b1' D1(t) x(t) + b2' D2(t) x(t) + u(t)

In all three the intercept shifts at the two break dates through a1 and a2, so a
break in the level is always present. In model(c) that is the only break:
the constant shifts at each break date while the slope b stays the same across
the whole sample. model(ct) keeps the same level shifts and adds a linear
trend g t; the trend itself does not break (there is no trend-by-dummy term), so
the constant breaks but the trend does not. In model(rs) both the constant
and the slopes break: the interaction terms D1(t) x(t) and D2(t) x(t) move the
slope from b0 to b0 + b1 after the first break and to b0 + b1 + b2 after the
second, so each regime has its own intercept and its own slope. model(rs)
is the default and is the model in Hatemi-J (2008, equation 2).

These tests extend the single-regime-shift cointegration tests of Gregory and
Hansen (1996) to two unknown regime shifts. The three deterministic
specifications follow Gregory and Hansen's model taxonomy (level shift; level
shift with trend; regime shift), here with two breaks rather than one, while
Hatemi-J (2008) develops the regime-shift case in equation (2) and supplies new
critical values for the ADF*, Zt*, and Za* statistics.

Cointegration is assessed through the residual u(t) of the chosen regression.
ADF* is the t-statistic on the lagged residual in the augmented Dickey-Fuller
regression of that residual; Zt* and Za* are the Phillips (1987) statistics
built from the bias-corrected first-order autoregressive coefficient of the
residual, using a long-run-variance kernel. Each statistic is taken as the
minimum over all admissible break pairs, the null hypothesis is no
cointegration, and more negative values provide evidence against it. The
reported critical values are those of Hatemi-J (2008) and cover up to four
regressors.

This command is a Stata port of Hatemi-J's coint2b routine of Gauss.

Options

model(string) sets the deterministic specification. model(rs) (regime
shift, the model in Hatemi-J 2008, eq. 2) is the default. model(c) uses a
constant with level shifts only; model(ct) adds a linear trend.

choice(#) sets how the ADF lag length is chosen. choice(2) (AIC) is
the default, matching Hatemi-J's coint2b routine of Gauss. choice(1)
uses kmax lags; choice(3) uses BIC; choice(4) uses the downward-t
procedure with a 1.96 cutoff on the longest lag. choice(5) uses the same
downward t-stat procedure but with the 1.645 cutoff.

choice(6) selects the lag by a general-to-specific Breusch-Godfrey (BG)
procedure that targets the elimination of serial correlation: starting from kmax
and stepping down, the ADF auxiliary regression is fitted at each lag and its
residuals are tested for autocorrelation with the BG LM test of orders 1 to
bglags(); the shortest lag whose residuals are still free of
autocorrelation (smallest BG p-value at or above 0.05) is kept. If
autocorrelation cannot be removed even at kmax, kmax is used and a warning is
printed. The header reports the smallest BG p-value over orders 1 to bglags()
at the selected lag; a value at or above 0.05 indicates no autocorrelation. The
rationale is that the information-criterion rules choice(2) (AIC) and
choice(3) (BIC) minimize a fit-versus-parsimony criterion and give no
guarantee that the residuals of the selected model are free of serial
correlation, while choice(4) and choice(5) only check the
significance of the longest lag. The asymptotic distributions underlying ADF*,
Zt* and Za* require approximately white-noise residuals, so a lag chosen by AIC
or BIC can leave residual autocorrelation that distorts the test.
choice(6) instead selects the lag so that the Breusch-Godfrey test detects
no residual autocorrelation, making the white-noise requirement an explicit
selection target rather than a hoped-for by-product.

kmax(#) sets the maximum number of lagged terms in the ADF
auxiliary regression. Default is 12.

kernel(string) chooses the long-run variance estimator used for the Zt* and
Za* statistics. kernel(qs) (quadratic spectral) is the default for the
coint2b method and reproduces Hatemi-J's code; kernel(iid) applies no
autocovariance correction; kernel(bartlett) uses the Bartlett kernel.

bwl(#) fixes the bandwidth (the number of autocovariance lags) used by the
bartlett and qs kernels. When it is not specified the default
bandwidth is the Andrews (1991) automatic value for the coint2b method
(1.3221(a*T)^(1/5) for qs, 1.1447(a*T)^(1/3) for bartlett, with a from
an AR(1) fit to the residuals) and round(4*(T/100)^(2/9)) for the tspdlib
method. Leaving bwl() unset therefore reproduces each reference exactly. The
option has no effect under kernel(iid).

trimming(#) sets the trimming rate t for the break search: the break dates
are searched over [t*T, (1-t)*T] and the two breaks must be at least t*T apart,
following Gregory and Hansen (1996). The default is 0.15, which reproduces
Hatemi-J's settings. The value must lie strictly between 0 and 0.5. This option
applies to both the default and the tspdlib method.

bglags(#) sets the highest Breusch-Godfrey order tested by choice(6)
(BG orders 1 to bglags() are all tested and the smallest p-value is used).
It has no effect under other lag-selection rules. When it is not specified the
default follows the data frequency of tsset, as in the kapetanios
command: 8 for quarterly, 24 for monthly, 52 for weekly, 100 for daily, and 2
otherwise (including annual series), capped at kmax*5.

reg(string) re-runs the cointegrating regression with regress so the
OLS estimates can be inspected directly. string selects the break used:
adf, zt, or za; a bare reg uses zt. It adds the dummies
D1, D2 and, under model rs, the interactions D1_var, D2_var to the
data, overwriting any like-named variables. The regress command line is
printed.

regnewey(string) is reg() with Newey-West HAC standard errors
(newey) and truncation lag floor(4*(T/100)^(2/9)). The coefficients match
OLS; only the standard errors differ. The newey command line is printed.

tspdlib computes the statistics with the TSPDLIB coint_hatemiJ routine
instead of coint2b.

Examples

Example 1. Replication of the output of Hatemi-J's coint2b routine of Gauss.

The data are the nelsonplosser series, regressing money on the bond yield:

. use "https://eruygurakademi.com/datasets/hatemij/nelsonplosser.dta", clear

. tsset year

. hatemij m bnd, model(rs) choice(2) kmax(12)

To reproduce the same numbers in Hatemi-J's GAUSS code, create a folder
named gaussexample on the C: drive (so its path is C:\gaussexample),
download the three files below into that folder, then open the .gss file in GAUSS
and run it (Run / F5). The GAUSS code is by Abdulnasser Hatemi-J (author of the
Hatemi-J test), available at
https://github.com/aptech/gauss-hatemij:

cItest2b.src

nelsonplosser.dta

hatemij_gauss_check.gss

Both produce the same results (model rs = regime shift, choice 2 = AIC,
kmax = 12):

Statistic Value

ADF* -6.490912 (lag 5; breaks 0.2676 / 0.4648)

Zt* -4.389183 (breaks 0.2394 / 0.4648)

Za* -34.188667 (breaks 0.2394 / 0.4648)

hatemij reproduces Hatemi-J's GAUSS output.

Example 2. A four-regressor regime-shift model: Error Case

. use "https://eruygurakademi.com/datasets/hatemij/nelsonplosser.dta", clear

. tsset year

. hatemij m bnd vel lrgnp lun, model(rs) choice(2) kmax(12)

To reproduce the same numbers in Hatemi-J's GAUSS code, create a folder
named gaussexample on the C: drive (so its path is C:\gaussexample),
download the three files below into that folder, then open the .gss file in GAUSS
and run it (Run / F5). The GAUSS code is by Abdulnasser Hatemi-J (author of the
Hatemi-J test), available at
https://github.com/aptech/gauss-hatemij:

cItest2b.src

nelsonplosser.dta

hatemij_gauss_error.gss

The same test is implemented in the TSPDLIB library as coint_hatemiJ.
To run it, create a folder named gaussexample on the C: drive (so its path
is C:\gaussexample), download the two files below into that folder, then open the
.gss file in GAUSS and run it (Run / F5). This requires the TSPDLIB library to be
installed in GAUSS (library tspdlib;). The call uses model 3 (regime shift,
corresponding to model(rs) here), the Akaike criterion, twelve lags, and a
trimming rate of 0.15 (the 15 percent trimming of Hatemi-J's code; the TSPDLIB
default is 0.10):

nelsonplosser.dta

hatemij_gauss_tspdlib.gss

With the same 15 percent trimming, all three implementations search the
same grid and locate the same minimizing break pair, at observations 30 and 39.
Hatemi-J's GAUSS code (coint2b) and the TSPDLIB routine (coint_hatemiJ)
both report the second break at observation 30, a value taken from the grid
column index rather than from the underlying break date; the two reported breaks
then coincide. For coint2b this makes the final cointegrating regression
singular (GAUSS returns a "matrix not positive definite" error). hatemij
reports the true second break at observation 39, so its cointegrating regression
is well defined and the command completes, returning ADF* = -7.985408,
Zt* = -8.051679, and Za* = -64.056591. The Zt* and Za* break locations agree
across the three implementations. The small differences in the Zt* and Za* values
come from the sample-size normalization each routine uses: coint2b forms
Za* as n times (rho-hat - 1) while TSPDLIB uses (n-1), and the two use n-1 versus
n-2 degrees of freedom in the Zt* variance. For this design the long-run-variance
correction is zero in both (the automatic bandwidth falls below one, so no
autocovariance terms enter), so the iid versus quadratic-spectral choice does not
contribute here.

Example 3. The tspdlib option

. use "https://eruygurakademi.com/datasets/hatemij/nelsonplosser.dta", clear

. tsset year

. hatemij m bnd vel lun, model(rs) choice(2) kmax(12) tspdlib

To reproduce the same numbers with the TSPDLIB library, create a folder
named gaussexample on the C: drive (so its path is C:\gaussexample),
download the two files below into that folder, then open the .gss file in GAUSS
and run it (Run / F5). This requires the TSPDLIB library to be installed in GAUSS
(library tspdlib;). The call uses model 3 (regime shift, corresponding to
model(rs) here), the Akaike criterion, twelve lags, the iid long-run
variance, and a trimming rate of 0.15 (the 15 percent trimming of Hatemi-J's
code; the TSPDLIB default is 0.10):

nelsonplosser.dta

hatemij_gauss_tspdlib3.gss

Both produce the same statistics and break dates:

Statistic Value

ADF* -7.178643 (lag 6; breaks 1916 / 1931)

Zt* -7.046005 (breaks 1942 / 1945)

Za* -58.417297 (breaks 1942 / 1945)

hatemij prints ADF* = -7.178643, Zt* = -7.046005, and
Za* = -58.417297. The TSPDLIB routine prints these same statistics rounded to
three decimals (-7.179, -7.046, and -58.417) and reports the identical break
dates, so the tspdlib option reproduces TSPDLIB's coint_hatemiJ.

Example 4. kernel and bwl options

. hatemij m bnd vel lrgnp, kernel(iid)

. hatemij m bnd vel lrgnp, kernel(bartlett) bwl(8)

. hatemij m bnd vel lrgnp, trimming(0.20)

Leaving kernel(), bwl(), and trimming() at their defaults
reproduces Hatemi-J's coint2b. The same options apply under tspdlib,
where they reproduce TSPDLIB's coint_hatemiJ with the corresponding
varm and bandwidth; trimming() applies to both methods.

Example 5. kernel and bwl options together with tspdlib option

. use "https://eruygurakademi.com/datasets/hatemij/nelsonplosser.dta", clear

. tsset year

. hatemij m bnd vel lun, model(rs) choice(2) kmax(12) tspdlib kernel(bartlett) bwl(8)

To reproduce the same numbers with the TSPDLIB library, create a folder
named gaussexample on the C: drive (so its path is C:\gaussexample),
download the two files below into that folder, then open the .gss file in GAUSS
and run it (Run / F5). This requires the TSPDLIB library to be installed in GAUSS
(library tspdlib;). The call uses model 3 (regime shift, corresponding to
model(rs) here), the Akaike criterion, twelve lags, the Bartlett long-run
variance (varm = 2), a bandwidth of 8, and a trimming rate of 0.15:

nelsonplosser.dta

hatemij_gauss_tspdlib_bart.gss

Both produce the same statistics and break dates:

Statistic Value

ADF* -7.178643 (lag 6; breaks 1916 / 1931)

Zt* -7.043270 (breaks 1942 / 1945)

Za* -58.260282 (breaks 1942 / 1945)

hatemij prints ADF* = -7.178643, Zt* = -7.043270, and
Za* = -58.260282. The TSPDLIB routine prints these same statistics rounded to
three decimals (-7.179, -7.043, and -58.260) and reports the identical break
dates. The ADF* statistic and its break dates are unchanged from the iid case
because the kernel affects only the Phillips statistics.

Example 6. The quadratic spectral kernel

. use "https://eruygurakademi.com/datasets/hatemij/nelsonplosser.dta", clear

. tsset year

. hatemij m bnd vel lun, model(rs) choice(2) kmax(12) tspdlib kernel(qs) bwl(8)

To reproduce the same numbers with the TSPDLIB library, create a folder
named gaussexample on the C: drive (so its path is C:\gaussexample),
download the two files below into that folder, then open the .gss file in GAUSS
and run it (Run / F5). This requires the TSPDLIB library to be installed in GAUSS
(library tspdlib;). The call uses model 3 (regime shift, corresponding to
model(rs) here), the Akaike criterion, twelve lags, the quadratic spectral
long-run variance (varm = 3), a bandwidth of 8, and a trimming rate of
0.15:

nelsonplosser.dta

hatemij_gauss_tspdlib_qs.gss

Both produce the same statistics and break dates:

Statistic Value

ADF* -7.178643 (lag 6; breaks 1916 / 1931)

Zt* -7.043502 (breaks 1942 / 1945)

Za* -58.273653 (breaks 1942 / 1945)

To full precision the TSPDLIB values are -7.043502018541 and
-58.273652989658; hatemij agrees with them to about ten decimal places. The
break dates coincide and ADF* is unchanged from the iid case.

Example 7. choice(5) option

. webuse lutkepohl2

. tsset qtr

. hatemij ln_inv ln_inc, model(rs) choice(5) kmax(8) kernel(iid)

hatemij reports:

Statistic Value

ADF* -5.516255 (lag 0; breaks 1966q3 / 1974q1, obs 27 / 57)

Zt* -6.047909 (breaks 1966q3 / 1973q4, obs 27 / 56)

Za* -53.127385 (breaks 1966q3 / 1973q4, obs 27 / 56)

ADF* and Za* match the TSPDLIB-style n-normalized results to the
precision those implementations report (-5.516 and -53.127).

A consistency check for Zt*. The Zt* and Za* statistics are both
functions of the same estimated first-order autoregressive coefficient at the
selected break, so they must agree in sign and order of magnitude. A large
negative Za* (here -53.13) cannot occur together with a near-zero Zt* at the same
break; a pair like that (for instance Za* = -53.13 reported with Zt* = -0.53)
signals an error in the Zt* computation rather than weak evidence against the
null. hatemij returns Zt* = -6.047909, consistent with Za* = -53.127385.

Example 8. choice(6) option: selecting the ADF model without autocorrelation

. webuse lutkepohl2

. tsset qtr

. hatemij ln_inv ln_inc, model(rs) choice(6) kmax(8) kernel(iid)

The header reports the lag rule and the Breusch-Godfrey result at the
selected lag, and hatemij returns:

Lag selection Breusch-Godfrey GTS (bglags = 8) (kmax = 8, selected lag = 0)

Breusch-Godfrey min p = 0.1905 (orders 1-8): no residual autocorrelation

Statistic Value

ADF* -5.548741 (lag 0; breaks 1965q4 / 1973q4, obs 24 / 56)

Zt* -6.047909 (breaks 1966q3 / 1973q4, obs 27 / 56)

Za* -53.127385 (breaks 1966q3 / 1973q4, obs 27 / 56)

The reported min p is the smallest Breusch-Godfrey p-value over orders 1 to
bglags() at the selected lag, evaluated at the ADF* break; here it is 0.19,
so the residuals show no autocorrelation. If autocorrelation persisted at every
lag up to kmax, the largest lag would be used and a warning would be printed
alongside the (significant) min p. The value is also stored in r(bgp).

Example 9. AIC Lag selection

. webuse lutkepohl2

. hatemij ln_inv ln_inc, model(rs) choice(2) kmax(12) tspdlib kernel(iid)

. hatemicoint ln_inv ln_inc, maxlags(12) lagselection(aic)

To reproduce the GAUSS results, create a folder named gaussexample on
the C: drive (path C:\gaussexample), download the two files below into it, and run
the .gss file in GAUSS (Run / F5) with the TSPDLIB library installed
(library tspdlib;); the call uses model 3 (regime shift), the default
bandwidth, the Akaike criterion (ic = 1), twelve lags, the iid long-run
variance (varm = 1), and a trimming rate of 0.15:

lutkepohl2.dta

hatemij_gauss_tspdlib5.gss

hatemij reports ADF* = -6.240, Zt* = -6.015 and Za* = -52.550.

tspdlib (GAUSS routine) reports ADF* = -6.240, Zt* = -6.015 and Za* = -52.550.

hatemicoint reports ADF* = -6.240, Zt* = -0.527 and Za* = -53.127.

Example 10. The quadratic spectral kernel with bandwidth 8

. webuse lutkepohl2

. hatemij ln_inv ln_inc, model(rs) choice(4) kmax(8) tspdlib kernel(qs) bwl(8)

. hatemicoint ln_inv ln_inc, kernel(qs) bwl(8)

To reproduce the GAUSS results, create a folder named gaussexample on
the C: drive (path C:\gaussexample), download the two files below into it, and run
the .gss file in GAUSS (Run / F5) with the TSPDLIB library installed
(library tspdlib;); the call uses model 3 (regime shift), bandwidth 8, the
t-stat lag rule (ic = 3), eight lags, the quadratic spectral long-run
variance (varm = 3), and a trimming rate of 0.15:

lutkepohl2.dta

hatemij_gauss_tspdlib4.gss

hatemij reports ADF* = -5.516, Zt* = -6.207 and Za* = -58.690.

tspdlib (GAUSS routine) reports ADF* = -5.516, Zt* = -6.207 and Za* = -58.690.

hatemicoint reports ADF* = -5.516, Zt* = -0.525 and Za* = -53.123.

Stored results

hatemij stores the following in r():

Scalars

r(adf) ADF* statistic (Modified ADF)

r(lag_adf) AR lag chosen at ADF*

r(tb1_adf) first break fraction (ADF*)

r(tb2_adf) second break fraction (ADF*)

r(obs1_adf) first break observation number (ADF*)

r(obs2_adf) second break observation number (ADF*)

r(date1_adf) first break time-variable value (ADF*)

r(date2_adf) second break time-variable value (ADF*)

r(zt) Zt* statistic (Modified Phillips)

r(tb1_zt) first break fraction (Zt*)

r(tb2_zt) second break fraction (Zt*)

r(obs1_zt) first break observation number (Zt*)

r(obs2_zt) second break observation number (Zt*)

r(date1_zt) first break time-variable value (Zt*)

r(date2_zt) second break time-variable value (Zt*)

r(za) Za* statistic (Modified Phillips)

r(tb1_za) first break fraction (Za*)

r(tb2_za) second break fraction (Za*)

r(obs1_za) first break observation number (Za*)

r(obs2_za) second break observation number (Za*)

r(date1_za) first break time-variable value (Za*)

r(date2_za) second break time-variable value (Za*)

r(N) number of observations

r(m) number of independent variables

r(bgp) minimum Breusch-Godfrey p-value at the selected lag (choice(6) only)

Matrices

r(ptable) parameter table (coef, std. err., t) at the Zt* break

r(cv_adfzt) 1%, 5%, 10% critical values for ADF*/Zt* (m <= 4)

r(cv_za) 1%, 5%, 10% critical values for Za* (m <= 4)

References

Andrews, D. W. K. 1991. Heteroskedasticity and autocorrelation consistent
covariance matrix estimation. Econometrica 59: 817-858.

Gregory, A. W., and B. E. Hansen. 1996. Residual-based tests for cointegration
in models with regime shifts. Journal of Econometrics 70: 99-126.

Hatemi-J, A. 2008. Tests for cointegration with two unknown regime shifts with
an application to financial market integration. Empirical Economics 35: 497-505.

Phillips, P. C. B. 1987. Time series regression with a unit root. Econometrica
55: 277-301.

Author

H. Ozan Eruygur

AHBV University, Ankara, Turkiye.

Department of Economics

https://www.ozaneruygur.com

eruygur@gmail.com

Eruygur Academy and Consulting (Eruygur Akademi ve Danismanlik), Ankara, Turkiye.

https://www.eruygurakademi.com

eruygurakademi@gmail.com

This command is a faithful Stata/Mata port of Hatemi-J's GAUSS code (coint2b).

hatemij v0.9.0 - June 2026

Please cite as:

Eruygur, H. O. 2026. hatemij: Hatemi-J (2008) cointegration test with two unknown regime shifts.
Stata package version 0.9.0. Available from: https://www.eruygurakademi.com