//import delimited "/Users/shanchao/OneDrive - University of California, Davis/KnowledgeStock/GammaModel/TwoSeries_testdata.csv", case(upper) 

tsset YEAR
gen LOG_MFP = log(MFP)
gen LOG_STOCK = log(STOCK)

// test stationary for the series and first difference
dfgls LOG_STOCK
matrix temp = r(cvalues)
gen DFGLS_LOG_STOCK_LAG = r(sclag) // number of SC lags
gen DFGLS_LOG_STOCK     = temp[10 - r(sclag) + 1, 2] //column 2 DF-GLS statistics
gen DFGLS_LOG_STOCK_1   = temp[10 - r(sclag) + 1, 3] //1% critical value
gen DFGLS_LOG_STOCK_5   = temp[10 - r(sclag) + 1, 4] //5% critical value
gen DFGLS_LOG_STOCK_10  = temp[10 - r(sclag) + 1, 5] //10% critical value

dfgls D.LOG_STOCK
matrix temp = r(cvalues)
gen DFGLS_D_LOG_STOCK_LAG = r(sclag) // number of SC lags
gen DFGLS_D_LOG_STOCK     = temp[10 - r(sclag) + 1, 2] //column 5 DF-GLS statistics
gen DFGLS_D_LOG_STOCK_1   = temp[10 - r(sclag) + 1, 3] //1% critical value
gen DFGLS_D_LOG_STOCK_5   = temp[10 - r(sclag) + 1, 4] //5% critical value
gen DFGLS_D_LOG_STOCK_10  = temp[10 - r(sclag) + 1, 5] //10% critical value

// testing cointegration PP
quietly reg LOG_MFP LOG_STOCK
quietly predict temp_resi, resid
pperron temp_resi, trend // no restriction
matrix temp = r(cvalues)
gen PPTEST_LAG  = r(lags)
gen PPTEST_T    = temp[2,1] //tau-statistics
gen PPTEST_T_1  = temp[2,2] //1% critical value 
gen PPTEST_T_5  = temp[2,3] //5% critical value
gen PPTEST_T_10 = temp[2,4] //10% critical value

drop temp_resi

// Johansen test
// lag order selection by AIC, HQIC, and SBIC
varsoc LOG_MFP LOG_STOCK, maxlag(10)
matrix temp = r(stats)
mata 
mata clear 
M = st_matrix("temp")
i = 0; w = 0; ii = 0; ww = 0; iii=0; www=0;
minindex(M[,7],1,i,w) // AIC
minindex(M[,8],1,ii,ww) // HQIC
minindex(M[,9],1,iii,www) // SBIC
res = (i - 1, ii - 1, iii-1)
st_matrix("J", res)
end

gen AIC  = J[1,1] // AIC lag
gen HQIC = J[1,2] // HQIC lag
gen SBIC = J[1,3] // SBIC lag

// AIC results
local AIC = J[1,1]
vecranknew LOG_MFP LOG_STOCK, lag(`AIC') trend(constant) // model without trend in cointegration equation 
matrix temp = e(trace)
gen AIC_TRACE_0_M3 = temp[1,1]
gen AIC_TRACE_1_M3 = temp[1,2]
gen AIC_TRACE_CV_0_M3 = critval0
gen AIC_TRACE_CV_1_M3 = critval1

local AIC = J[1,1]
vecranknew LOG_MFP LOG_STOCK, lag(`AIC') trend(rtrend) // model with trend in cointegration equation 
matrix temp = e(trace)
gen AIC_TRACE_0_M4 = temp[1,1]
gen AIC_TRACE_1_M4 = temp[1,2]
gen AIC_TRACE_CV_0_M4 = critval0
gen AIC_TRACE_CV_1_M4 = critval1


// HQIC results
local HQIC = J[1,2]
vecranknew LOG_MFP LOG_STOCK, lag(`HQIC') trend(constant) // model without trend in cointegration equation 
matrix temp = e(trace)
gen HQIC_TRACE_0_M3 = temp[1,1]
gen HQIC_TRACE_1_M3 = temp[1,2]
gen HQIC_TRACE_CV_0_M3 = critval0
gen HQIC_TRACE_CV_1_M3 = critval1

local HQIC = J[1,2]
vecranknew LOG_MFP LOG_STOCK, lag(`HQIC') trend(rtrend) // model with trend in cointegration equation 
matrix temp = e(trace)
gen HQIC_TRACE_0_M4 = temp[1,1]
gen HQIC_TRACE_1_M4 = temp[1,2]
gen HQIC_TRACE_CV_0_M4 = critval0
gen HQIC_TRACE_CV_1_M4 = critval1

// SBIC results
local SBIC = J[1,3]
vecranknew LOG_MFP LOG_STOCK, lag(`SBIC') trend(constant) // model without trend in cointegration equation 
matrix temp = e(trace)
gen SBIC_TRACE_0_M3 = temp[1,1]
gen SBIC_TRACE_1_M3 = temp[1,2]
gen SBIC_TRACE_CV_0_M3 = critval0
gen SBIC_TRACE_CV_1_M3 = critval1

local SBIC = J[1,3]
vecranknew LOG_MFP LOG_STOCK, lag(`SBIC') trend(rtrend) // model with trend in cointegration equation 
matrix temp = e(trace)
gen SBIC_TRACE_0_M4 = temp[1,1]
gen SBIC_TRACE_1_M4 = temp[1,2]
gen SBIC_TRACE_CV_0_M4 = critval0
gen SBIC_TRACE_CV_1_M4 = critval1

gen id = _n
keep if id == 1
drop YEAR INDEX STOCK MFP TREND LOG_MFP LOG_STOCK id