'''
[first posted version November 2009]
Code used for the co-breaking analysis in:
Schreiber, Sven (2009)
Unemployment and Productivity, Slowdowns and Speed-Ups; Evidence Using Common Shifts
Published in The B.E. Journal of Macroeconomics, 2009, Vol. 9 : Iss. 1 (Topics), Article 39.
www.bepress.com/bejm/vol9/iss1/art39
DOI: 10.2202/1935-1690.1818.
License: same as Numerical Python (NumPy), but please cite the source (the paper) if you use it.
Instructions:
To use this code you need Python (2.x, x>4 I think) with NumPy and SciPy installed.
You also need the module 'helpers.py' in the same directory as this file (you should be able to get it from http://econ.schreiberlin.de/software/helpers.py).
Finally, you need to put your data in a comma-separated-values (.csv) file, put that file in the same directory and specify the name below in the FILENAME variable.
The data file needs the following structure: the actual variables in the first columns, then the shift (step) dummies, then the impulse dummies. (You need to include the corresponding impulse dummies for the shift dates by hand, but you could also include additional impulses, which is the reason why the impulse dummy specification is not automated.)
If it doesn't work you can contact me at or the address provided on econ.schreiberlin.de (if that's different).
Caveat: there is still some clumsy hard-coding in the beginning, but it's clearly marked and you shouldn't have much problems substituting it.
'''
#### User input
# In my application the ordering of the variables was:
# ur_sa,prodgrowth,s74q4,s86q4,s96q2,i74q4,i86q4,i96q2
# (where the "s" variables are step dummies, and the "i" variables are impulse dummies)
FILENAME = 'yourdata.csv'
NUMOFVARS = 2
NUMOFBREAKS = 3 # must be >= NUMOFVARS
MAXLAG = 3
CBRANK = 1 # this is NUMOFVARS - (number of shift packs = rank of shift coeffs)
#### end of user input
from numpy.matlib import ( r_, c_, mat, ones, zeros, empty, rand, eye, log)
from numpy.linalg import eig, eigh, solve, lstsq, det
from scipy import stats as ST
from helpers import ( geneigsympos, addLags, null, getDeterministics,
vecm2varcoeffs, readcsv, vec, unvec, write_gretl_mat_xml )
### --- in the following lines there is hardcoding ---
data, orientation, varnames, obslabels, freq = readcsv(FILENAME)
endowlags = addLags(data[:,:2], MAXLAG)
# (this '2' should probably be substituted by NUMOFVARS)
impulseswlags = addLags(data[:,-3:], MAXLAG)
# (this '3' refers to the number of impulses, so number of cols in .csv file minus NUMOFVARS minus NUMOFBREAKS)
shifts = data[MAXLAG:, 2:5]
# (this '2': NUMOFVARS; this '5': NUMOFVARS + NUMOFBREAKS)
#regime_74 = 1 - shifts[:,0] # (benchmark regime)
# ok here it should really have been a loop to define the regime dummies
regime74_86 = shifts[:,0] - shifts[:,1]
regime86_96 = shifts[:,1] - shifts[:,2]
regime96_end = shifts[:,2]
#regimes = c_[regime74_86, regime86_96, regime96_end] # unrestricted version
regimes = c_[regime74_86, regime86_96] # restricted version
### --- end hardcoding ---
endo = endowlags[:,:NUMOFVARS] # endogenous variables
teff = endo.shape[0] # effective sample size
endolags_cnst_impul = c_[endowlags[:,NUMOFVARS:], ones((teff,1)), impulseswlags]
# concentrating:
R0 = endo - endolags_cnst_impul * lstsq(endolags_cnst_impul, endo)[0]
R1 = regimes - endolags_cnst_impul * lstsq(endolags_cnst_impul, regimes)[0]
# moment matrices
M00 = R0.T * R0 / teff
M01 = R0.T * R1 / teff
M11 = R1.T * R1 / teff
M = M01.T * solve(M00, M01)
# generalized eigenvalues, corresponding to det(lambda_j M11 - M)=0
# need the evals ascending first for tracestats
evals, evectors = geneigsympos(M, M11)
tracestats = -teff * log(1-evals).cumsum()
tracestats = tracestats[::-1] # reversed
evals = evals[::-1] # accordingly descending
evectors = evectors[:, ::-1]
print '-------- Co-breaking rank test --------------'
print 'rank -- eigenvalue -- trace stat -- p-value'
# (in scipy.stats the method sf stands for survival function = pvalue = 1-cdf)
for b in range(NUMOFVARS):
dof = (NUMOFVARS-b) * (NUMOFBREAKS-b)
print b, evals[b], tracestats[b], ST.chi2.sf(tracestats[b], dof)
print '---------------------------------------------'
# estimate the coefficients for pre-specified CBRANK
#print 'all eigenvectors'
#print evectors
beta = evectors[:,:NUMOFVARS-CBRANK]
# ok let's normalize this on the first break (second regime):
beta /= beta[0].A1
print 'The normalized shift coeff. matrix, beta:'
print beta
commonshifts = regimes * beta
restcoeffs, res, rk = lstsq(c_[endolags_cnst_impul, commonshifts], endo)[:3]
# The dim of this should be ...x NUMOFVARS:
# # 1stlag*NUMOFVARS, 2ndlag*NUMOFVARS, ..., MAXLAGthlag*NUMOFVARS
# # one constant term
# # NUMOFBREAKS * (MAXLAG+1) impulse dummies
# # (NUMOFVARS-CBRANK) shift packages
# So the constant term coefficient from the reduced-rank-UVAR, k0 aka kappa:
k0T = restcoeffs[MAXLAG*NUMOFVARS] # should be row vector now
print 'The conditional means k0T:', k0T
# the loadings of the common shift packages in the RR-UVAR, alpha aka xi:
alphaT = restcoeffs[-(NUMOFVARS-CBRANK):]
print 'The loadings in the UVAR, alphaT:', alphaT
# \Gamma(1) = I - \sum_{i=1}^g \Gamma_i to transform the estimates
Gamma1T = eye(NUMOFVARS)
for lag in range(MAXLAG):
Gamma1T -= restcoeffs[lag * NUMOFVARS : (lag+1) * NUMOFVARS]
print 'Eigenvalues of G1T:'
temp = eig(Gamma1T)
print temp[0]
# derive the coefficients of the UC form
G1TI = Gamma1T.I
c0T = k0T * G1TI
deltaT = alphaT * G1TI
deltaorthT = null(deltaT.T).T
print '--------------------------'
print 'The loadings in the UC form, deltaT:'
print deltaT
print 'Co-breaking vector normalized on the first variable:'
deltaorthT /= deltaorthT[:,0].A1 # .A1 is necessary to get broadcasting
print deltaorthT
print 'The "unconditional" means of the variables:'
# so endo * deltaorthT.T should be stationary!
# and the overall mean function:
print '1st regime (pre74):', c0T
for regime in range(regimes.shape[1]):
print 'regime ' + str(regime+2) + ' :', c0T + beta[regime] * deltaT
print '...and the mean of the co-breaking combination:'
print c0T * deltaorthT.T
print '-----------------------------------'
print '-----------------------------------'