From 1529d864693ad34b6e9980fecc7e4003106382c0 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Thu, 2 Jul 2020 15:40:48 +0800
Subject: [PATCH 01/28] Update accounting.py

Update some characteristics
Change the definition of ME
Add prcc_c into comp
Use Linear Regression to get ir
---
 pychars/accounting.py | 56 ++++++++++++++++++++++++++++++++++++-------
 1 file changed, 47 insertions(+), 9 deletions(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index 299fff9..ba4eb3f 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -5,7 +5,6 @@
 from dateutil.relativedelta import *
 from pandas.tseries.offsets import *
 import pickle as pkl
-
 ###################
 # Connect to WRDS #
 ###################
@@ -72,7 +71,7 @@ def ttm12(series, df):
                     f.ceq, f.scstkc, f.emp, f.csho, f.seq, f.txditc, f.pstkrv, f.pstkl, f.np, f.txdc, f.dpc, f.ajex,
                     
                     /*market*/
-                    abs(f.prcc_f) as prcc_f
+                    abs(f.prcc_f) as prcc_f, abs(f.prcc_c) as prcc_c
                     
                     from comp.funda as f
                     left join comp.company as c
@@ -210,7 +209,7 @@ def ttm12(series, df):
 ccm1 = pd.merge(comp, ccm, how='left', on=['gvkey'])
 
 # we can only get the accounting data after the firm public their report
-# for annual data, we ues 6 months lagged data
+# for annual data, we use 6 months lagged data
 ccm1['yearend'] = ccm1['datadate'] + YearEnd(0)
 ccm1['jdate'] = ccm1['yearend'] + MonthEnd(6)
 
@@ -218,6 +217,7 @@ def ttm12(series, df):
 ccm2 = ccm1[(ccm1['jdate'] >= ccm1['linkdt']) & (ccm1['jdate'] <= ccm1['linkenddt'])]
 
 # link comp and crsp
+# data_rawa only includes annul data because comp is annual. Use inner merge
 crsp2 = crsp2.rename(columns={'monthend': 'jdate'})
 data_rawa = pd.merge(crsp2, ccm2, how='inner', on=['permno', 'jdate'])
 
@@ -229,8 +229,8 @@ def ttm12(series, df):
 '''
 Note: me is CRSP market equity, mve_f is Compustat market equity. Please choose the me below.
 '''
-# data_rawa['me'] = data_rawa['me']/1000  # CRSP ME
-data_rawa['me'] = data_rawa['mve_f']  # Compustat ME
+data_rawa['me'] = data_rawa['me']/1000  # CRSP ME
+#data_rawa['me'] = data_rawa['mve_f']  # Compustat ME
 
 # update count after merging
 # data_rawa['count'] = data_rawa.groupby(['gvkey']).cumcount() + 1
@@ -289,10 +289,48 @@ def ttm12(series, df):
               np.nan]
 data_rawa['cfp_n'] = np.select(condlist, choicelist, default=data_rawa['ib']+data_rawa['dp'])
 
-# ep
+# ep, checked from Hou and change 'ME' from compustat to crsp
 data_rawa['ep'] = data_rawa['ib']/data_rawa['me']
 data_rawa['ep_n'] = data_rawa['ib']
 
+#ir
+'''
+First calculate r(t-5,t). Then rb(t-5,t) and use Bm to perform linear regression and get residue
+'''
+#r(t-5,t):sum ret from t-5 to t (which is calendar year t-6 to t-1)
+lag = pd.DataFrame()
+for i in range(1,6):
+        lag['ret%s' % i] = data_rawa.groupby(['permno'])['ret'].shift(i)
+        log_ps['ps%s'% i] = 
+data_rawa['ret5'] = lag['ret1']+lag['ret2']+lag['ret3']+lag['ret4']+lag['ret5']
+
+#bm_t-5 (bm of year t-5)
+data_rawa['bm5'] = data_rawa.groupby(['permno'])['bm'].shift(5)
+
+#rB (five year log book return)
+#Reference: jf_06 page8 by KENT DANIEL
+data_rawa['rB'] = data_rawa['bm'] - data_rawa['bm5'] + data_rawa['ret5']
+
+#Regression and get ir
+#First get unique datelist
+datelist = data_rawa['jdate'].unique()
+for date in datelist:
+    temp = data_rawa['jdate' == date]
+    n_row = temp.shape[0]
+    index = temp.index
+    X = pd.DataFrame()
+    X['bm5'] = temp['bm5']
+    X['rB'] = temp['rB']
+    X['intercept'] = 1
+    X = X[['intercept','rB','bm5']]
+    X = np.mat(X)
+    Y = np.mat(temp[['ret5']])
+    #These are residuals on one date
+    res = (np.identity(n_row) - X.dot(X.T.dot(X).I).dot(X.T)).dot(Y)
+    #put residuals back into data_rawa
+    data_rawa.loc[index,'ir'] = res
+
+
 # ni
 data_rawa['csho_l1'] = data_rawa.groupby(['permno'])['csho'].shift(1)
 data_rawa['ajex_l1'] = data_rawa.groupby(['permno'])['ajex'].shift(1)
@@ -524,8 +562,8 @@ def ttm12(series, df):
 '''
 Note: me is CRSP market equity, mveq_f is Compustat market equity. Please choose the me below.
 '''
-# data_rawq['me'] = data_rawq['me']/1000  # CRSP ME
-data_rawq['me'] = data_rawq['mveq_f']  # Compustat ME
+data_rawq['me'] = data_rawq['me']/1000  # CRSP ME
+#data_rawq['me'] = data_rawq['mveq_f']  # Compustat ME
 
 # update count after merging
 # data_rawq['count'] = data_rawq.groupby(['gvkey']).cumcount() + 1
@@ -594,7 +632,7 @@ def ttm12(series, df):
                           (ttm4('ibq', data_rawq)+ttm4('dpq', data_rawq))/data_rawq['me'])
 data_rawq['cfp_n'] = data_rawq['cfp']*data_rawq['me']
 
-# ep
+# ep, also checked and change 'ME' from compustat to crsp
 data_rawq['ep'] = ttm4('ibq', data_rawq)/data_rawq['me']
 data_rawq['ep_n'] = data_rawq['ep']*data_rawq['me']
 

From fc0e1db428ace6699405193ac7e0875769d82dc2 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Tue, 7 Jul 2020 11:01:17 +0800
Subject: [PATCH 02/28] Update accounting.py

---
 pychars/accounting.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index ba4eb3f..559f3ba 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -71,7 +71,7 @@ def ttm12(series, df):
                     f.ceq, f.scstkc, f.emp, f.csho, f.seq, f.txditc, f.pstkrv, f.pstkl, f.np, f.txdc, f.dpc, f.ajex,
                     
                     /*market*/
-                    abs(f.prcc_f) as prcc_f, abs(f.prcc_c) as prcc_c
+                    abs(f.prcc_f) as prcc_f, abs(f.prcc_c) as prcc_c, f.dvc, f.prstkc, f.sstk, f.fopt, f.wcap, f.oancf
                     
                     from comp.funda as f
                     left join comp.company as c
@@ -339,7 +339,7 @@ def ttm12(series, df):
                            np.log(data_rawa['csho']*data_rawa['ajex']).replace(-np.inf, 0)-
                            np.log(data_rawa['csho_l1']*data_rawa['ajex_l1']).replace(-np.inf, 0))
 
-# op
+# op: the formula seems different from Hou Page 74?
 data_rawa['cogs0'] = np.where(data_rawa['cogs'].isnull(), 0, data_rawa['cogs'])
 data_rawa['xint0'] = np.where(data_rawa['xint'].isnull(), 0, data_rawa['xint'])
 data_rawa['xsga0'] = np.where(data_rawa['xsga'].isnull(), 0, data_rawa['xsga'])
@@ -349,6 +349,18 @@ def ttm12(series, df):
 data_rawa['op'] = np.select(condlist, choicelist,
                           default=(data_rawa['revt'] - data_rawa['cogs0'] - data_rawa['xsga0'] - data_rawa['xint0'])/data_rawa['be'])
 
+#nop
+data_rawa['net_p'] = data_rawa['dvc'] + data_rawa['prstkc'] + 2*data_rawa['pstkrv'] - data_rawa['sstk']
+data_rawa['nop'] = data_rawa['net_p'] / data_rawa['me']
+data_rawa['nop'] = np.where(data_rawa['nop']<=0, nan, data_rawa['nop'] )
+
+#ocp
+data_rawa['ocy'] = np.where(data_rawa['jdate'] < '1988-06-30', data_rawa['fopt'] - data_rawa['wcap'], data_rawa['fopt'] - data_rawa['oancf'])
+data_rawa['ocp'] = data_rawa['ocy'] / data_rawa['me']
+data_rawa['ocp'] = np.where(data_rawa['ocp']<=0, nan, data_rawa['ocp'] )
+
+
+
 # rsup
 data_rawa['sale_l1'] = data_rawa.groupby(['permno'])['sale'].shift(1)
 data_rawa['rsup'] = (data_rawa['sale']-data_rawa['sale_l1'])/data_rawa['me']
@@ -363,7 +375,7 @@ def ttm12(series, df):
 # lev
 data_rawa['lev'] = data_rawa['lt']/data_rawa['me']
 
-# sp
+# sp, checked
 data_rawa['sp'] = data_rawa['sale']/data_rawa['me']
 data_rawa['sp_n'] = data_rawa['sale']
 

From 85f1783e89be6a4d99c4ae467a304d939ef4c374 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Tue, 14 Jul 2020 20:58:03 +0800
Subject: [PATCH 03/28] Update accounting.py

---
 pychars/accounting.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index be899a8..af5f9b1 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -270,7 +270,7 @@ def ttm12(series, df):
               np.nan]
 data_rawa['cfp_n'] = np.select(condlist, choicelist, default=data_rawa['ib']+data_rawa['dp'])
 
-# ep, checked from Hou and change 'ME' from compustat to crsp
+# ep, checked from Hou and change 'ME' from compustat to crsp，checked
 data_rawa['ep'] = data_rawa['ib']/data_rawa['me']
 data_rawa['ep_n'] = data_rawa['ib']
 
@@ -340,8 +340,6 @@ def ttm12(series, df):
 data_rawa['ocp'] = data_rawa['ocy'] / data_rawa['me']
 data_rawa['ocp'] = np.where(data_rawa['ocp']<=0, nan, data_rawa['ocp'] )
 
-
-
 # rsup
 data_rawa['sale_l1'] = data_rawa.groupby(['permno'])['sale'].shift(1)
 data_rawa['rsup'] = (data_rawa['sale']-data_rawa['sale_l1'])/data_rawa['me']

From c1a615f964ccd4bbba8c74e32b9b90526f75593d Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Mon, 20 Jul 2020 17:58:20 +0800
Subject: [PATCH 04/28] Update accounting.py

---
 pychars/accounting.py | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index af5f9b1..b45005a 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -282,7 +282,7 @@ def ttm12(series, df):
 lag = pd.DataFrame()
 for i in range(1,6):
         lag['ret%s' % i] = data_rawa.groupby(['permno'])['ret'].shift(i)
-        log_ps['ps%s'% i] = 
+
 data_rawa['ret5'] = lag['ret1']+lag['ret2']+lag['ret3']+lag['ret4']+lag['ret5']
 
 #bm_t-5 (bm of year t-5)
@@ -481,6 +481,16 @@ def ttm12(series, df):
 # dy
 data_rawa['dy'] = data_rawa['dvt']/data_rawa['me']
 
+#aci
+data_rawa['capx_s'] = data_rawa['capx']/data_rawa['sale']
+data_rawa['capx_3'] = data_rawa.groupby('permno')['capx_s'].shift(1)+data_rawa.groupby('permno')['capx_s'].shift(2)+data_rawa.groupby('permno')['capx_s'].shift(3)
+data_rawa['capx_3'] = data_rawa['capx_3']/3
+data_rawa['aci'] = data_rawa['capx_s']/data_rawa['capx_3']
+
+#cei
+data_rawa['me_5'] = data_rawa.groupby('permno')['me'].shift(5)
+data_rawa['cei'] = np.log(data_rawa['me']/data_rawa['me_5']) - data_rawa['ret5']
+
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',
                      'sic', 'acc', 'agr', 'bm', 'cfp', 'ep', 'ni', 'op', 'rsup', 'cash', 'chcsho',

From 7276644edff26088103c253f707b6b1a84b30e19 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Tue, 21 Jul 2020 11:48:29 +0800
Subject: [PATCH 05/28] Update chars

---
 pychars/accounting.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index b45005a..3e95b6e 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -491,6 +491,12 @@ def ttm12(series, df):
 data_rawa['me_5'] = data_rawa.groupby('permno')['me'].shift(5)
 data_rawa['cei'] = np.log(data_rawa['me']/data_rawa['me_5']) - data_rawa['ret5']
 
+#dwc
+data_rawa['dwc'] = (data_rawa['act'] - data_rawa['che']) - (data_rawa['lct'] - data_rawa['dlc'])
+
+#I/A
+data_rawa['ia'] = (data_rawa['at']/data_rawa['at_l1'])-1
+
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',
                      'sic', 'acc', 'agr', 'bm', 'cfp', 'ep', 'ni', 'op', 'rsup', 'cash', 'chcsho',

From 477e37a0f897b992078ebcf87b14d1f6a4581221 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Wed, 22 Jul 2020 11:15:44 +0800
Subject: [PATCH 06/28] Update accounting.py

---
 pychars/accounting.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index 3e95b6e..b02feba 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -778,6 +778,12 @@ def ttm12(series, df):
 # ato
 data_rawq['ato'] = data_rawq['saleq']/data_rawq['noa_l4']
 
+#Iaq
+data_rawq['atqlag'] = ttm4('atq',data_rawq)
+data_rawq['iaq'] = (data_rawq['atq']/data_rawq['atqlag'])-1
+
+
+
 # Quarterly Accounting Variables
 chars_q = data_rawq[['gvkey', 'permno', 'datadate', 'jdate', 'sic', 'exchcd', 'shrcd', 'acc', 'bm', 'cfp',
                      'ep', 'agr', 'ni', 'op', 'cash', 'chcsho', 'rd', 'cashdebt', 'pctacc', 'gma', 'lev',
@@ -881,4 +887,6 @@ def mom(start, end, df):
     pkl.dump(chars_a, f)
 
 with open('chars_q.pkl', 'wb') as f:
-    pkl.dump(chars_q, f)
\ No newline at end of file
+    pkl.dump(chars_q, f)
+
+

From 8b3eead27b608883b8112ef3011ce64cfc92fe5a Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Thu, 23 Jul 2020 11:22:16 +0800
Subject: [PATCH 07/28] Update another 4 investment chars

---
 pychars/accounting.py | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index b02feba..0076f72 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -65,7 +65,7 @@ def ttm12(series, df):
                     
                     /*liabilities*/
                     f.lct, f.dlc, f.dltt, f.lt, f.dm, f.dcvt, f.cshrc, 
-                    f.dcpstk, f.pstk, f.ap, f.lco, f.lo, f.drc, f.drlt, f.txdi,
+                    f.dcpstk, f.pstk, f.ap, f.lco, f.lo, f.drc, f.drlt, f.txdi, f.dltis, f.dltr. f.dlcch
                     
                     /*equity and other*/
                     f.ceq, f.scstkc, f.emp, f.csho, f.seq, f.txditc, f.pstkrv, f.pstkl, f.np, f.txdc, f.dpc, f.ajex,
@@ -497,6 +497,21 @@ def ttm12(series, df):
 #I/A
 data_rawa['ia'] = (data_rawa['at']/data_rawa['at_l1'])-1
 
+#Ig
+data_rawa['capx_l1'] = data_rawa.groupby('permno')['capx'].shift(1)
+data_rawa['ig'] = data_rawa['capx']/data_rawa['capx_l1']
+
+#2Ig
+data_rawa['capx_l2'] = data_rawa.groupby('permno')['capx'].shift(2)
+data_rawa['ig'] = data_rawa['capx']/data_rawa['capx_l2']
+
+#Ivc
+data_rawa['atAvg'] = (data_rawa['at']+data_rawa['at_l1'])/2
+data_rawa['ivc'] = data_rawa['invt'] / data_rawa['atAvg']
+
+#Ndf
+data_rawa['ndf'] = data_rawa['dltis'] - data_rawa['dltr'] + data_rawa['dlcch'] 
+
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',
                      'sic', 'acc', 'agr', 'bm', 'cfp', 'ep', 'ni', 'op', 'rsup', 'cash', 'chcsho',

From db73ca9469188549730272b8eb3313f786e29fe1 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Thu, 23 Jul 2020 18:12:59 +0800
Subject: [PATCH 08/28] update 2 chars

---
 pychars/accounting.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index 0076f72..dd5961f 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -419,7 +419,7 @@ def ttm12(series, df):
 # alm
 # data_rawa['alm'] = data_rawa['ala']/(data_rawa['at']+data_rawa['prcc_f']*data_rawa['csho']-data_rawa['ceq'])
 
-# noa
+# noa,checked
 data_rawa['noa'] = ((data_rawa['at']-data_rawa['che']-data_rawa['ivao'].fillna(0))-
                   (data_rawa['at']-data_rawa['dlc'].fillna(0)-data_rawa['dltt'].fillna(0)-data_rawa['mib'].fillna(0)
                    -data_rawa['pstk'].fillna(0)-data_rawa['ceq'])/data_rawa['at_l1'])
@@ -512,6 +512,12 @@ def ttm12(series, df):
 #Ndf
 data_rawa['ndf'] = data_rawa['dltis'] - data_rawa['dltr'] + data_rawa['dlcch'] 
 
+#nsi
+data_rawa['sps'] = data_rawa['csho'] * data_rawa['ajex']
+data_rawa['sps_l1'] = data_rawa.groupby('permno')['sps'].shift(1)
+data_rawa['nsi'] = np.log(data_rawa['sps']/data_rawa['sps_l1'])
+
+
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',
                      'sic', 'acc', 'agr', 'bm', 'cfp', 'ep', 'ni', 'op', 'rsup', 'cash', 'chcsho',

From 62adbb9bbeb0d1e095961e0baa8449e5a66b460b Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Fri, 24 Jul 2020 19:45:47 +0800
Subject: [PATCH 09/28] update until ta except pta

---
 pychars/accounting.py | 24 +++++++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index dd5961f..9b4a288 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -58,7 +58,7 @@ def ttm12(series, df):
                     f.ebit, f.nopi, f.spi, f.pi, f.txp, f.ni, f.txfed, f.txfo, f.txt, f.xint,
                     
                     /*CF statement and others*/
-                    f.capx, f.oancf, f.dvt, f.ob, f.gdwlia, f.gdwlip, f.gwo, f.mib, f.oiadp, f.ivao,
+                    f.capx, f.oancf, f.dvt, f.ob, f.gdwlia, f.gdwlip, f.gwo, f.mib, f.oiadp, f.ivao, f.ivst
                     
                     /*assets*/
                     f.rect, f.act, f.che, f.ppegt, f.invt, f.at, f.aco, f.intan, f.ao, f.ppent, f.gdwl, f.fatb, f.fatl,
@@ -493,6 +493,7 @@ def ttm12(series, df):
 
 #dwc
 data_rawa['dwc'] = (data_rawa['act'] - data_rawa['che']) - (data_rawa['lct'] - data_rawa['dlc'])
+#data_rawa['dwc'] = data_rawa['dwc']/data_rawa['at_l1']
 
 #I/A
 data_rawa['ia'] = (data_rawa['at']/data_rawa['at_l1'])-1
@@ -517,6 +518,27 @@ def ttm12(series, df):
 data_rawa['sps_l1'] = data_rawa.groupby('permno')['sps'].shift(1)
 data_rawa['nsi'] = np.log(data_rawa['sps']/data_rawa['sps_l1'])
 
+#oa
+data_rawa['txp'] = np.where(data_rawa['txp'].isnull(), 0, data_rawa['txp'])
+data_rawa['oa'] = (data_rawa['act'] - data_rawa['che']) - (data-rawa['lct'] - data_rawa['dlc'] - data_rawa['txp']) - data_rawa['dp']
+
+#Poa
+data_rawa['poa'] = data_rawa['oa']/data_rawa['ni']
+
+#dNco
+data_rawa['lct'] = np.where(data_rawa['lct'].isnull(), 0, data_rawa['lct'])
+data_rawa['dltt'] = np.where(data_rawa['dltt'].isnull(), 0, data_rawa['dltt'])
+data_rawa['ivao'] = np.where(data_rawa['ivao'].isnull(), 0, data_rawa['ivao'])
+data_rawa['ivst'] = np.where(data_rawa['ivst'].isnull(), 0, data_rawa['ivst'])
+data_rawa['pstk'] = np.where(data_rawa['pstk'].isnull(), 0, data_rawa['pstk'])
+data_rawa['dnco'] = (data_rawa['at'] - data_rawa['ivao']) - (data_rawa['lt'] - data_rawa['dltt'])
+
+#dFin
+data_rawa['dfin'] = (data_rawa['ivst'] + data_rawa['ivao']) - (data_rawa['dltt'] + data_rawa['dlc'] + data_rawa['pstk'])
+
+#Ta
+data_rawa['ta'] = data_rawa['dwc'] + data_rawa['dnco'] + data_rawa['dfin']
+
 
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',

From 95601f8d2ca0fdd8aa042c273975a96cc518dc6f Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Mon, 27 Jul 2020 18:26:32 +0800
Subject: [PATCH 10/28] fixed bugs

---
 pychars/accounting.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index 9b4a288..ac9d5c9 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -58,14 +58,14 @@ def ttm12(series, df):
                     f.ebit, f.nopi, f.spi, f.pi, f.txp, f.ni, f.txfed, f.txfo, f.txt, f.xint,
                     
                     /*CF statement and others*/
-                    f.capx, f.oancf, f.dvt, f.ob, f.gdwlia, f.gdwlip, f.gwo, f.mib, f.oiadp, f.ivao, f.ivst
+                    f.capx, f.oancf, f.dvt, f.ob, f.gdwlia, f.gdwlip, f.gwo, f.mib, f.oiadp, f.ivao, f.ivst,
                     
                     /*assets*/
                     f.rect, f.act, f.che, f.ppegt, f.invt, f.at, f.aco, f.intan, f.ao, f.ppent, f.gdwl, f.fatb, f.fatl,
                     
                     /*liabilities*/
                     f.lct, f.dlc, f.dltt, f.lt, f.dm, f.dcvt, f.cshrc, 
-                    f.dcpstk, f.pstk, f.ap, f.lco, f.lo, f.drc, f.drlt, f.txdi, f.dltis, f.dltr. f.dlcch
+                    f.dcpstk, f.pstk, f.ap, f.lco, f.lo, f.drc, f.drlt, f.txdi, f.dltis, f.dltr, f.dlcch,
                     
                     /*equity and other*/
                     f.ceq, f.scstkc, f.emp, f.csho, f.seq, f.txditc, f.pstkrv, f.pstkl, f.np, f.txdc, f.dpc, f.ajex,
@@ -333,12 +333,12 @@ def ttm12(series, df):
 #nop
 data_rawa['net_p'] = data_rawa['dvc'] + data_rawa['prstkc'] + 2*data_rawa['pstkrv'] - data_rawa['sstk']
 data_rawa['nop'] = data_rawa['net_p'] / data_rawa['me']
-data_rawa['nop'] = np.where(data_rawa['nop']<=0, nan, data_rawa['nop'] )
+data_rawa['nop'] = np.where(data_rawa['nop']<=0, np.nan, data_rawa['nop'] )
 
 #ocp
 data_rawa['ocy'] = np.where(data_rawa['jdate'] < '1988-06-30', data_rawa['fopt'] - data_rawa['wcap'], data_rawa['fopt'] - data_rawa['oancf'])
 data_rawa['ocp'] = data_rawa['ocy'] / data_rawa['me']
-data_rawa['ocp'] = np.where(data_rawa['ocp']<=0, nan, data_rawa['ocp'] )
+data_rawa['ocp'] = np.where(data_rawa['ocp']<=0, np.nan, data_rawa['ocp'] )
 
 # rsup
 data_rawa['sale_l1'] = data_rawa.groupby(['permno'])['sale'].shift(1)
@@ -483,9 +483,9 @@ def ttm12(series, df):
 
 #aci
 data_rawa['capx_s'] = data_rawa['capx']/data_rawa['sale']
-data_rawa['capx_3'] = data_rawa.groupby('permno')['capx_s'].shift(1)+data_rawa.groupby('permno')['capx_s'].shift(2)+data_rawa.groupby('permno')['capx_s'].shift(3)
-data_rawa['capx_3'] = data_rawa['capx_3']/3
-data_rawa['aci'] = data_rawa['capx_s']/data_rawa['capx_3']
+data_rawa['capx3'] = data_rawa.groupby(['permno'])['capx_s'].shift(1) + data_rawa.groupby(['permno'])['capx_s'].shift(2) + data_rawa.groupby(['permno'])['capx_s'].shift(3)
+data_rawa['capx3'] = data_rawa['capx3']/3
+data_rawa['aci'] = data_rawa['capx_s']/data_rawa['capx3']
 
 #cei
 data_rawa['me_5'] = data_rawa.groupby('permno')['me'].shift(5)
@@ -504,7 +504,7 @@ def ttm12(series, df):
 
 #2Ig
 data_rawa['capx_l2'] = data_rawa.groupby('permno')['capx'].shift(2)
-data_rawa['ig'] = data_rawa['capx']/data_rawa['capx_l2']
+data_rawa['2ig'] = data_rawa['capx']/data_rawa['capx_l2']
 
 #Ivc
 data_rawa['atAvg'] = (data_rawa['at']+data_rawa['at_l1'])/2
@@ -520,7 +520,7 @@ def ttm12(series, df):
 
 #oa
 data_rawa['txp'] = np.where(data_rawa['txp'].isnull(), 0, data_rawa['txp'])
-data_rawa['oa'] = (data_rawa['act'] - data_rawa['che']) - (data-rawa['lct'] - data_rawa['dlc'] - data_rawa['txp']) - data_rawa['dp']
+data_rawa['oa'] = (data_rawa['act'] - data_rawa['che']) - (data_rawa['lct'] - data_rawa['dlc'] - data_rawa['txp']) - data_rawa['dp']
 
 #Poa
 data_rawa['poa'] = data_rawa['oa']/data_rawa['ni']

From 1b435ab3c264cc6953a0b0e14b94df48aa2276a5 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Tue, 28 Jul 2020 17:57:06 +0800
Subject: [PATCH 11/28] adm and almq

---
 pychars/accounting.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index ac9d5c9..fbb689e 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -364,7 +364,7 @@ def ttm12(series, df):
 # rdm
 data_rawa['rdm'] = data_rawa['xrd']/data_rawa['me']
 
-# adm hxz adm
+# adm hxz adm, checked
 data_rawa['adm'] = data_rawa['xad']/data_rawa['me']
 
 # gma
@@ -540,6 +540,7 @@ def ttm12(series, df):
 data_rawa['ta'] = data_rawa['dwc'] + data_rawa['dnco'] + data_rawa['dfin']
 
 
+
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',
                      'sic', 'acc', 'agr', 'bm', 'cfp', 'ep', 'ni', 'op', 'rsup', 'cash', 'chcsho',
@@ -825,7 +826,12 @@ def ttm12(series, df):
 data_rawq['atqlag'] = ttm4('atq',data_rawq)
 data_rawq['iaq'] = (data_rawq['atq']/data_rawq['atqlag'])-1
 
-
+#Almq
+data_rawq['intanq'] = np.where(data_rawq['intanq'].isnull(), 0, data_rawq['intanq'])
+data_rawq['qal'] = data_rawq['cheq'] + 0.75*(data_rawq['actq']-data_rawq['cheq']) + 0.5*(data_rawq['atq'] - data_rawq['actq'] - data_rawq['intanq'])
+data_rawq['mveqa'] = data_rawq['atq'] + data_rawq['mveq_f'] - data_rawq['ceqq']
+data_rawq['mveqa_1'] = data_rawq.groupby(['permno'])['mveqa'].shift(1)
+data_rawq['almq'] = data_rawq['qal']/data_rawq['mveqa_1']
 
 # Quarterly Accounting Variables
 chars_q = data_rawq[['gvkey', 'permno', 'datadate', 'jdate', 'sic', 'exchcd', 'shrcd', 'acc', 'bm', 'cfp',

From 28c8a3bd01fe2053e3933e0fbe209423f5d567c2 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Thu, 30 Jul 2020 10:05:49 +0800
Subject: [PATCH 12/28] ol,etr

---
 pychars/accounting.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index fbb689e..34e9bb9 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -68,7 +68,7 @@ def ttm12(series, df):
                     f.dcpstk, f.pstk, f.ap, f.lco, f.lo, f.drc, f.drlt, f.txdi, f.dltis, f.dltr, f.dlcch,
                     
                     /*equity and other*/
-                    f.ceq, f.scstkc, f.emp, f.csho, f.seq, f.txditc, f.pstkrv, f.pstkl, f.np, f.txdc, f.dpc, f.ajex,
+                    f.ceq, f.scstkc, f.emp, f.csho, f.seq, f.txditc, f.pstkrv, f.pstkl, f.np, f.txdc, f.dpc, f.ajex, f.epspx,
                     
                     /*market*/
                     abs(f.prcc_f) as prcc_f, abs(f.prcc_c) as prcc_c, f.dvc, f.prstkc, f.sstk, f.fopt, f.wcap, f.oancf
@@ -539,7 +539,16 @@ def ttm12(series, df):
 #Ta
 data_rawa['ta'] = data_rawa['dwc'] + data_rawa['dnco'] + data_rawa['dfin']
 
+#Ol
+data_rawa['ol'] = (data_rawa['cogs'] + data_rawa['xsga'])/data_rawa['at']
 
+#etr
+data_rawa['txtpi'] = data_rawa['txt'] / data_rawa['pi']
+data_rawa['txtpi_l1'] = data_rawa.groupby('permno')['txtpi'].shift(1)
+data_rawa['txtpi_l2'] = data_rawa.groupby('permno')['txtpi'].shift(2)
+data_rawa['txtpi_l3'] = data_rawa.groupby('permno')['txtpi'].shift(3)
+data_rawa['deps'] = data_rawa['epspx']/(data_rawa['ajex'] * data_rawa['prcc_f'])
+data_rawa['etr'] = (data_rawa['txtpi'] - (data_rawa['txtpi_l1'] + data_rawa['txtpi_l2'] + data_rawa['txtpi_l3'])/3) * data_rawa['deps']
 
 # Annual Accounting Variables
 chars_a = data_rawa[['cusip', 'ncusip', 'gvkey', 'permno', 'exchcd', 'shrcd', 'datadate', 'jdate', 'count',

From 2f0d589c0bb665dd4494dd350fbc55427e0c8dca Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Fri, 31 Jul 2020 10:41:48 +0800
Subject: [PATCH 13/28] olq and seasonality

---
 pychars/accounting.py | 82 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 82 insertions(+)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index 34e9bb9..6745160 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -842,6 +842,10 @@ def ttm12(series, df):
 data_rawq['mveqa_1'] = data_rawq.groupby(['permno'])['mveqa'].shift(1)
 data_rawq['almq'] = data_rawq['qal']/data_rawq['mveqa_1']
 
+#Olq
+data_rawa['olq'] = (data_rawa['cogsq'] + data_rawa['xsgaq'])/data_rawa['atq']
+
+
 # Quarterly Accounting Variables
 chars_q = data_rawq[['gvkey', 'permno', 'datadate', 'jdate', 'sic', 'exchcd', 'shrcd', 'acc', 'bm', 'cfp',
                      'ep', 'agr', 'ni', 'op', 'cash', 'chcsho', 'rd', 'cashdebt', 'pctacc', 'gma', 'lev',
@@ -862,6 +866,84 @@ def ttm12(series, df):
 crsp_mom['jdate'] = pd.to_datetime(crsp_mom['date']) + MonthEnd(0)
 crsp_mom = crsp_mom.dropna()
 
+#Seasonality
+
+#Rla
+crsp_mom['rla'] = crsp_mom.groupby['permno']['ret'].shift(12)
+
+#Rln
+lag = pd.DataFrame()
+result = 0
+for i in range(1, 12):
+    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
+    result = result + lag['mom%s' % i]
+crsp_mom['rln'] = result/11
+
+#R[2,5]a
+#R[2,5]n
+lag = pd.DataFrame()
+result = 0
+for i in range(13,61):
+    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
+    if i not in [24,36,48,60]:
+        result = result + lag['mom%s' % i]
+
+crsp_mom['r25a'] = (lag['mom24']+lag['mom36']+lag['mom48']+lag['mom60'])/4
+crsp_mom['r25n'] = result/44
+
+#R[6,10]a
+#R[6,10]n
+lag = pd.DataFrame()
+result = 0
+for i in range(61,121):
+    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
+    if i not in [72,84,96,108,120]:
+        result = result + lag['mom%s' % i]
+
+crsp_mom['r610a'] = (lag['mom72']+lag['mom84']+lag['mom96']+lag['mom108']+lag['mom120'])/5
+crsp_mom['r610n'] = result/55
+
+#R[11,15]a
+lag = pd.DataFrame()
+result = 0
+for i in [132,144,156,168,180]:
+    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
+    result = result + lag['mom%s' % i]
+crsp_mom['r1115a'] = result/5
+
+#R[16,20]a
+lag = pd.DataFrame()
+result = 0
+for i in [192,204,216,228,240]:
+    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
+    result = result + lag['mom%s' % i]
+crsp_mom['r1620a'] = result/5
+
+#R[2,5]n
+lag = pd.DataFrame()
+result = 0
+for i in [24,36,48,60]:
+    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
+    result = result + lag['mom%s' % i]
+crsp_mom['rln'] = result/4
+
+
+def mom(start, end, df):
+    """
+
+    :param start: Order of starting lag
+    :param end: Order of ending lag
+    :param df: Dataframe
+    :return: Momentum factor
+    """
+    lag = pd.DataFrame()
+    result = 1
+    for i in range(start, end):
+        lag['mom%s' % i] = df.groupby(['permno'])['ret'].shift(i)
+        result = result * (1+lag['mom%s' % i])
+    result = result - 1
+    return result
+
 # add delisting return
 dlret = conn.raw_sql("""
                      select permno, dlret, dlstdt 

From 7ced2be5618a674750be705a58b5ee6968f303e3 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Tue, 4 Aug 2020 09:06:17 +0800
Subject: [PATCH 14/28] rdm,rdmq,rds

---
 pychars/accounting.py | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/pychars/accounting.py b/pychars/accounting.py
index 6745160..d252b85 100644
--- a/pychars/accounting.py
+++ b/pychars/accounting.py
@@ -747,6 +747,9 @@ def ttm12(series, df):
 # rdm
 data_rawq['rdm'] = data_rawq['xrdq4']/data_rawq['me']
 
+# rds
+data_rawq['rds'] = data_rawq['xrdq4']/data_rawq['saleq']
+
 # sgr
 data_rawq['saleq4'] = ttm4('saleq', data_rawq)
 data_rawq['saleq4'] = np.where(data_rawq['saleq4'].isnull(), data_rawq['saley'], data_rawq['saleq4'])
@@ -919,13 +922,6 @@ def ttm12(series, df):
     result = result + lag['mom%s' % i]
 crsp_mom['r1620a'] = result/5
 
-#R[2,5]n
-lag = pd.DataFrame()
-result = 0
-for i in [24,36,48,60]:
-    lag['mom%s' % i] = crsp_mom.groupby(['permno'])['ret'].shift(i)
-    result = result + lag['mom%s' % i]
-crsp_mom['rln'] = result/4
 
 
 def mom(start, end, df):

From 6d58f56194e6f5421382ce75c916b441e490f3c3 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Mon, 10 Aug 2020 15:55:25 +0800
Subject: [PATCH 15/28] Create dtv.py

---
 pychars/dtv.py | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 74 insertions(+)
 create mode 100644 pychars/dtv.py

diff --git a/pychars/dtv.py b/pychars/dtv.py
new file mode 100644
index 0000000..f678b2e
--- /dev/null
+++ b/pychars/dtv.py
@@ -0,0 +1,74 @@
+# Fama & French 3 factors residual variance
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block
+crsp = conn.raw_sql("""
+                      select a.permno, a.date, a.ret, a.prc, a.vol
+                      from crsp.dsf as a
+                      where a.date > '01/01/1959'
+                      """)
+
+# sort variables by permno and date
+crsp = crsp.sort_values(by=['permno', 'date'])
+
+# change variable format to int
+crsp['permno'] = crsp['permno'].astype(int)
+
+# Line up date to be end of month
+crsp['date'] = pd.to_datetime(crsp['date'])
+
+# find the closest trading day to the end of the month
+crsp['monthend'] = crsp['date'] + MonthEnd(0)
+crsp['date_diff'] = crsp['monthend'] - crsp['date']
+date_temp = crsp.groupby(['permno', 'monthend'])['date_diff'].min()
+date_temp = pd.DataFrame(date_temp)  # convert Series to DataFrame
+date_temp.reset_index(inplace=True)
+date_temp.rename(columns={'date_diff': 'min_diff'}, inplace=True)
+crsp = pd.merge(crsp, date_temp, how='left', on=['permno', 'monthend'])
+crsp['sig'] = np.where(crsp['date_diff'] == crsp['min_diff'], 1, np.nan)
+
+# label every date of month end
+crsp['month_count'] = crsp[crsp['sig'] == 1].groupby(['permno']).cumcount()
+
+# mark the number of each month to each day of this month
+crsp['month_count'] = crsp.groupby(['permno'])['month_count'].fillna(method='bfill')
+
+
+######################
+#    Calculate dtv   #
+######################
+
+def main(start, end, step):
+    crsp['dtv'] = crsp['prc'] * crsp['vol']
+    return crsp
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+with open('rvar_ff3.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
\ No newline at end of file

From fb510d08a8424d172dd4f7f9525d44f47530de07 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Fri, 28 Aug 2020 22:28:52 +0800
Subject: [PATCH 16/28] create capm_mom

---
 pychars/{dtv.py => hxz_dtv.py} |   0
 pychars/rvar_capm_mom.py       | 168 +++++++++++++++++++++++++++++++++
 2 files changed, 168 insertions(+)
 rename pychars/{dtv.py => hxz_dtv.py} (100%)
 create mode 100644 pychars/rvar_capm_mom.py

diff --git a/pychars/dtv.py b/pychars/hxz_dtv.py
similarity index 100%
rename from pychars/dtv.py
rename to pychars/hxz_dtv.py
diff --git a/pychars/rvar_capm_mom.py b/pychars/rvar_capm_mom.py
new file mode 100644
index 0000000..926d21c
--- /dev/null
+++ b/pychars/rvar_capm_mom.py
@@ -0,0 +1,168 @@
+# CAPM residual momentum
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block. We use crsp.msf and ff.factors_monthly
+crsp = conn.raw_sql("""
+                      select a.permno, a.date, a.ret, (a.ret - b.rf) as exret, b.mktrf
+                      from crsp.msf as a
+                      left join ff.factors_monthly as b
+                      on a.date=b.date
+                      where a.date >= '01/01/1959'
+                      """)
+
+# sort variables by permno and date
+crsp = crsp.sort_values(by=['permno', 'date'])
+
+# change variable format to int
+crsp['permno'] = crsp['permno'].astype(int)
+
+# create a firm list. Use reset_index to create a count to number companies without duplicates
+df_firm = crsp.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+#Extract number of data points for each permno
+crsp['monthcount'] = crsp.groupby('permno').cumcount()
+monthnum = crsp.groupby('permno')['monthcount'].tail(1)
+monthnum = monthnum.astype(int)
+
+######################
+# Calculate residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    """
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    #for every permno, we have count as its number of obervations, prog as its number in firmlist
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        #this is to demonstrate how many companies have already been calculated, and convert it into percentage
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        #Actually there are count+1 months
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            # Now we have temp as sixty month data
+            temp = df[(df['permno'] == firm) & (i -59 <= df['month_count']) & (df['month_count'] <= i)]
+            # if observations in last 3 months are less 21, we drop the rvar of this month
+            if temp['permno'].count() < 21:
+                pass
+            else:
+                rolling_window = temp['permno'].count()
+                index = temp.tail(1).index
+                X = pd.DataFrame()
+                X[['mktrf']] = temp[['mktrf']]
+                X['intercept'] = 1
+                X = X[['intercept', 'mktrf']]
+                X = np.mat(X)
+                Y = np.mat(temp[['exret']])
+                res = (np.identity(rolling_window) - X.dot(X.T.dot(X).I).dot(X.T)).dot(Y)
+                df.loc[index, 'res'] = res
+    return df
+
+
+    def sub_df(start, end, step):
+    """
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+def main(start, end, step):
+    """
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_res_var, (df['crsp%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+# process dataframe
+crsp = crsp.dropna(subset=['res'])  # drop NA due to rolling
+crsp = crsp.rename(columns={'res': 'capm_1-1'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'capm_1-1']]
+
+def mom(start, end, df):
+    """
+    :param start: Order of starting lag
+    :param end: Order of ending lag
+    :param df: Dataframe
+    :return: Momentum factor
+    """
+    lag = pd.DataFrame()
+    result = 1
+    for i in range(start, end):
+        lag['mom%s' % i] = df.groupby(['permno'])['capm_1-1'].shift(i)
+        result = result * (1+lag['mom%s' % i])
+    result = result - 1
+    return result
+
+
+crsp['capm_2-12'] = mom(2,12,crsp)
+crsp['capm_13-60'] = mom(13,60,crsp)
+
+with open('rvar_capm_mom.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
+
+

From 6af48732d64878e651d6f29376b94a9ec51ad69a Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Sat, 29 Aug 2020 11:04:05 +0800
Subject: [PATCH 17/28] Update rvar_capm.py

---
 pychars/rvar_capm.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pychars/rvar_capm.py b/pychars/rvar_capm.py
index 08ec24c..847c2e1 100644
--- a/pychars/rvar_capm.py
+++ b/pychars/rvar_capm.py
@@ -86,7 +86,7 @@ def get_res_var(df, firm_list):
             # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
             temp = df[(df['permno'] == firm) & (i - 2 <= df['month_count']) & (df['month_count'] <= i)]
             # if observations in last 3 months are less 21, we drop the rvar of this month
-            if temp['permno'].count() < 21:
+            if temp['permno'].count() < 20:
                 pass
             else:
                 rolling_window = temp['permno'].count()
@@ -129,7 +129,6 @@ def sub_df(start, end, step):
 
 def main(start, end, step):
     """
-
     :param start: the quantile to start cutting, usually it should be 0
     :param end: the quantile to end cutting, usually it should be 1
     :param step: quantile step

From f09353c5db2370165b35544deeb12c0f6d6d8dcb Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Sat, 29 Aug 2020 15:25:38 +0800
Subject: [PATCH 18/28] Update rvar_capm_mom.py

Fix bugs
---
 pychars/rvar_capm_mom.py | 64 ++++++++++++++++------------------------
 1 file changed, 26 insertions(+), 38 deletions(-)

diff --git a/pychars/rvar_capm_mom.py b/pychars/rvar_capm_mom.py
index 926d21c..d418152 100644
--- a/pychars/rvar_capm_mom.py
+++ b/pychars/rvar_capm_mom.py
@@ -23,19 +23,24 @@
 conn = wrds.Connection()
 
 # CRSP Block. We use crsp.msf and ff.factors_monthly
-crsp = conn.raw_sql("""
-                      select a.permno, a.date, a.ret, (a.ret - b.rf) as exret, b.mktrf
-                      from crsp.msf as a
-                      left join ff.factors_monthly as b
-                      on a.date=b.date
+cr = conn.raw_sql("""
+                      select a.permno, a.date, a.ret from crsp.msf as a
                       where a.date >= '01/01/1959'
                       """)
 
-# sort variables by permno and date
-crsp = crsp.sort_values(by=['permno', 'date'])
+ff = conn.raw_sql("""
+                      select b.rf, b.mktrf, b.date from ff.factors_monthly as b
+                      where b.date >= '01/01/1959'
+                      """)
 
-# change variable format to int
+ff['date'] = pd.to_datetime(ff['date']) + MonthEnd(0)
+cr['date'] = pd.to_datetime(cr['date']) + MonthEnd(0)
+
+crsp = pd.merge(cr,ff,how = 'left', on = ['date'])
+crsp['exret'] = crsp['ret'] - crsp['rf']
+crsp = crsp.sort_values(by=['permno', 'date'])
 crsp['permno'] = crsp['permno'].astype(int)
+crsp = crsp[['permno','date','ret','exret','mktrf']]
 
 # create a firm list. Use reset_index to create a count to number companies without duplicates
 df_firm = crsp.drop_duplicates(['permno'])
@@ -44,11 +49,11 @@
 df_firm = df_firm.reset_index(drop=True)
 df_firm = df_firm.reset_index()
 df_firm = df_firm.rename(columns={'index': 'count'})
-
+print(df_firm)
 #Extract number of data points for each permno
-crsp['monthcount'] = crsp.groupby('permno').cumcount()
-monthnum = crsp.groupby('permno')['monthcount'].tail(1)
-monthnum = monthnum.astype(int)
+crsp['month_count'] = crsp.groupby('permno').cumcount()
+month_num = crsp.groupby('permno')['month_count'].tail(1)
+month_num = month_num.astype(int)
 
 ######################
 # Calculate residual #
@@ -56,11 +61,6 @@
 
 
 def get_res_var(df, firm_list):
-    """
-    :param df: stock dataframe
-    :param firm_list: list of firms matching stock dataframe
-    :return: dataframe with variance of residual
-    """
     #for every permno, we have count as its number of obervations, prog as its number in firmlist
     for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
         prog = prog + 1
@@ -72,11 +72,12 @@ def get_res_var(df, firm_list):
             # Now we have temp as sixty month data
             temp = df[(df['permno'] == firm) & (i -59 <= df['month_count']) & (df['month_count'] <= i)]
             # if observations in last 3 months are less 21, we drop the rvar of this month
-            if temp['permno'].count() < 21:
+            if temp['permno'].count() < 20:
                 pass
             else:
                 rolling_window = temp['permno'].count()
                 index = temp.tail(1).index
+                print(index)
                 X = pd.DataFrame()
                 X[['mktrf']] = temp[['mktrf']]
                 X['intercept'] = 1
@@ -84,17 +85,12 @@ def get_res_var(df, firm_list):
                 X = np.mat(X)
                 Y = np.mat(temp[['exret']])
                 res = (np.identity(rolling_window) - X.dot(X.T.dot(X).I).dot(X.T)).dot(Y)
-                df.loc[index, 'res'] = res
+                #print(df)
+                df.loc[index, 'res'] = res[-1]
     return df
 
 
-    def sub_df(start, end, step):
-    """
-    :param start: the quantile to start cutting, usually it should be 0
-    :param end: the quantile to end cutting, usually it should be 1
-    :param step: quantile step
-    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
-    """
+def sub_df(start, end, step):
     # we use dict to store different sub dataframe
     temp = {}
     for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
@@ -111,12 +107,6 @@ def sub_df(start, end, step):
     return temp
 
 def main(start, end, step):
-    """
-    :param start: the quantile to start cutting, usually it should be 0
-    :param end: the quantile to end cutting, usually it should be 1
-    :param step: quantile step
-    :return: a dataframe with calculated variance of residual
-    """
     df = sub_df(start, end, step)
     pool = mp.Pool()
     p_dict = {}
@@ -138,18 +128,13 @@ def main(start, end, step):
     crsp = main(0, 1, 0.05)
 
 # process dataframe
+print(crsp)
 crsp = crsp.dropna(subset=['res'])  # drop NA due to rolling
 crsp = crsp.rename(columns={'res': 'capm_1-1'})
 crsp = crsp.reset_index(drop=True)
 crsp = crsp[['permno', 'date', 'capm_1-1']]
 
 def mom(start, end, df):
-    """
-    :param start: Order of starting lag
-    :param end: Order of ending lag
-    :param df: Dataframe
-    :return: Momentum factor
-    """
     lag = pd.DataFrame()
     result = 1
     for i in range(start, end):
@@ -166,3 +151,6 @@ def mom(start, end, df):
     pkl.dump(crsp, f)
 
 
+
+
+

From 219f79ce0807b615ed2f1e5b0548af9a99ff98be Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Wed, 2 Sep 2020 21:57:35 +0800
Subject: [PATCH 19/28] Update rvar_capm_mom.py

---
 pychars/rvar_capm_mom.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pychars/rvar_capm_mom.py b/pychars/rvar_capm_mom.py
index d418152..a27a88e 100644
--- a/pychars/rvar_capm_mom.py
+++ b/pychars/rvar_capm_mom.py
@@ -49,7 +49,7 @@
 df_firm = df_firm.reset_index(drop=True)
 df_firm = df_firm.reset_index()
 df_firm = df_firm.rename(columns={'index': 'count'})
-print(df_firm)
+
 #Extract number of data points for each permno
 crsp['month_count'] = crsp.groupby('permno').cumcount()
 month_num = crsp.groupby('permno')['month_count'].tail(1)
@@ -77,7 +77,7 @@ def get_res_var(df, firm_list):
             else:
                 rolling_window = temp['permno'].count()
                 index = temp.tail(1).index
-                print(index)
+              
                 X = pd.DataFrame()
                 X[['mktrf']] = temp[['mktrf']]
                 X['intercept'] = 1
@@ -144,8 +144,8 @@ def mom(start, end, df):
     return result
 
 
-crsp['capm_2-12'] = mom(2,12,crsp)
-crsp['capm_13-60'] = mom(13,60,crsp)
+crsp['capm_2-12'] = mom(1,12,crsp)
+crsp['capm_13-60'] = mom(12,60,crsp)
 
 with open('rvar_capm_mom.pkl', 'wb') as f:
     pkl.dump(crsp, f)

From 2a0f1fefd2a5c967c6066925c6524e82d6c7b175 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Wed, 2 Sep 2020 22:03:30 +0800
Subject: [PATCH 20/28] Update rmom_capm

---
 pychars/{rvar_capm_mom.py => rmom_capm.py} | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)
 rename pychars/{rvar_capm_mom.py => rmom_capm.py} (96%)

diff --git a/pychars/rvar_capm_mom.py b/pychars/rmom_capm.py
similarity index 96%
rename from pychars/rvar_capm_mom.py
rename to pychars/rmom_capm.py
index a27a88e..06f9b00 100644
--- a/pychars/rvar_capm_mom.py
+++ b/pychars/rmom_capm.py
@@ -128,11 +128,10 @@ def main(start, end, step):
     crsp = main(0, 1, 0.05)
 
 # process dataframe
-print(crsp)
 crsp = crsp.dropna(subset=['res'])  # drop NA due to rolling
-crsp = crsp.rename(columns={'res': 'capm_1-1'})
+crsp = crsp.rename(columns={'res': 'rmom_capm_1m'})
 crsp = crsp.reset_index(drop=True)
-crsp = crsp[['permno', 'date', 'capm_1-1']]
+crsp = crsp[['permno', 'date', 'rmom_capm_1m']]
 
 def mom(start, end, df):
     lag = pd.DataFrame()
@@ -144,10 +143,10 @@ def mom(start, end, df):
     return result
 
 
-crsp['capm_2-12'] = mom(1,12,crsp)
-crsp['capm_13-60'] = mom(12,60,crsp)
+crsp['rmom_capm_12m'] = mom(1,12,crsp)
+crsp['rmom_capm_60m'] = mom(12,60,crsp)
 
-with open('rvar_capm_mom.pkl', 'wb') as f:
+with open('rmom_capm.pkl', 'wb') as f:
     pkl.dump(crsp, f)
 
 

From 6db51aba6948a01def5b83ff0a6b76042460b983 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Thu, 3 Sep 2020 01:03:04 +0800
Subject: [PATCH 21/28] Update rmom_capm.py

---
 pychars/rmom_capm.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pychars/rmom_capm.py b/pychars/rmom_capm.py
index 06f9b00..fcdc706 100644
--- a/pychars/rmom_capm.py
+++ b/pychars/rmom_capm.py
@@ -137,7 +137,7 @@ def mom(start, end, df):
     lag = pd.DataFrame()
     result = 1
     for i in range(start, end):
-        lag['mom%s' % i] = df.groupby(['permno'])['capm_1-1'].shift(i)
+        lag['mom%s' % i] = df.groupby(['permno'])['rmom_capm_1m'].shift(i)
         result = result * (1+lag['mom%s' % i])
     result = result - 1
     return result

From 8728451952704ad4f77fff7e999347a82d4f5cf8 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Tue, 8 Sep 2020 12:27:08 +0800
Subject: [PATCH 22/28] Create rmom_ff3.py

---
 pychars/rmom_ff3.py | 198 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 198 insertions(+)
 create mode 100644 pychars/rmom_ff3.py

diff --git a/pychars/rmom_ff3.py b/pychars/rmom_ff3.py
new file mode 100644
index 0000000..2fb3bc5
--- /dev/null
+++ b/pychars/rmom_ff3.py
@@ -0,0 +1,198 @@
+# Fama & French 3 factors residual variance
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block. We use crsp.msf and ff.factors_monthly
+cr = conn.raw_sql("""
+                      select a.permno, a.date, a.ret from crsp.msf as a
+                      where a.date >= '01/01/1959'
+                      """)
+
+ff = conn.raw_sql("""
+                      select b.rf, b.mktrf, b.date,b.smb, b.hml from ff.factors_monthly as b
+                      where b.date >= '01/01/1959'
+                      """)
+
+ff['date'] = pd.to_datetime(ff['date']) + MonthEnd(0)
+cr['date'] = pd.to_datetime(cr['date']) + MonthEnd(0)
+
+crsp = pd.merge(cr,ff,how = 'left', on = ['date'])
+crsp['exret'] = crsp['ret'] - crsp['rf']
+crsp = crsp.sort_values(by=['permno', 'date'])
+crsp['permno'] = crsp['permno'].astype(int)
+crsp = crsp[['permno','date','ret','exret','mktrf','smb','hml']]
+
+# create a firm list. Use reset_index to create a count to number companies without duplicates
+df_firm = crsp.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+#Extract number of data points for each permno
+crsp['month_count'] = crsp.groupby('permno').cumcount()
+month_num = crsp.groupby('permno')['month_count'].tail(1)
+month_num = month_num.astype(int)
+
+######################
+# Calculate the beta #
+######################
+# function that get multiple beta
+''''
+rolling_window = 60  # 60 trading days
+crsp['beta_mktrf'] = np.nan
+crsp['beta_smb'] = np.nan
+crsp['beta_hml'] = np.nan
+
+
+def get_beta(df):
+    """
+    The original idea of calculate beta is using formula (X'MX)^(-1)X'MY,
+    where M = I - 1(1'1)^{-1}1, I is a identity matrix.
+
+    """
+    temp = crsp.loc[df.index]  # extract the rolling sub dataframe from original dataframe
+    X = np.mat(temp[['mktrf', 'smb', 'hml']])
+    Y = np.mat(temp[['exret']])
+    ones = np.mat(np.ones(rolling_window)).T
+    M = np.identity(rolling_window) - ones.dot((ones.T.dot(ones)).I).dot(ones.T)
+    beta = (X.T.dot(M).dot(X)).I.dot((X.T.dot(M).dot(Y)))
+    crsp['beta_mktrf'].loc[df.index[-1:]] = beta[0]
+    crsp['beta_smb'].loc[df.index[-1:]] = beta[1]
+    crsp['beta_hml'].loc[df.index[-1:]] = beta[2]
+    return 0  # we do not need the rolling outcome since rolling cannot return different values in different columns
+
+
+# calculate beta through rolling window
+crsp_temp = crsp.groupby('permno').rolling(rolling_window).apply(get_beta, raw=False)
+'''
+
+######################
+# Calculate residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    """
+
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            temp = df[(df['permno'] == firm) & (i - 59 <= df['month_count']) & (df['month_count'] <= i)]
+            # if observations in last 3 months are less than 20, we drop the rvar of this month
+            if temp['permno'].count() < 20:
+                pass
+            else:
+                rolling_window = temp['permno'].count()
+                index = temp.tail(1).index
+                X = pd.DataFrame()
+                X[['mktrf', 'smb', 'hml']] = temp[['mktrf', 'smb', 'hml']]
+                X['intercept'] = 1
+                X = X[['intercept', 'mktrf', 'smb', 'hml']]
+                X = np.mat(X)
+                Y = np.mat(temp[['exret']])
+                res = (np.identity(rolling_window) - X.dot(X.T.dot(X).I).dot(X.T)).dot(Y)
+                df.loc[index, 'res'] = res[-1]
+    return df
+
+
+def sub_df(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+
+def main(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_res_var, (df['crsp%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+# process dataframe
+crsp = crsp.dropna(subset=['res'])  # drop NA due to rolling
+crsp = crsp.rename(columns={'res': 'rmom_ff3_1m'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'rmom_ff3_1m']]
+
+def mom(start, end, df):
+    lag = pd.DataFrame()
+    result = 1
+    for i in range(start, end):
+        lag['mom%s' % i] = df.groupby(['permno'])['rmom_ff3_1m'].shift(i)
+        result = result * (1+lag['mom%s' % i])
+    result = result - 1
+    return result
+
+
+crsp['rmom_ff3_12m'] = mom(1,12,crsp)
+crsp['rmom_ff3_60m'] = mom(12,60,crsp)
+
+with open('rmom_ff3.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
\ No newline at end of file

From 2029821fa49528ac59e70f7c03d0c6213625dcfe Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Mon, 21 Sep 2020 22:15:53 +0800
Subject: [PATCH 23/28] Create dtv.py

Calculate dtv from crsp.dsf
---
 pychars/dtv.py | 133 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 133 insertions(+)
 create mode 100644 pychars/dtv.py

diff --git a/pychars/dtv.py b/pychars/dtv.py
new file mode 100644
index 0000000..1f0e409
--- /dev/null
+++ b/pychars/dtv.py
@@ -0,0 +1,133 @@
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import pickle as pkl
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+#######################################################################################################################
+#                                                  CRSP Block                                                    #
+#######################################################################################################################
+
+# Create a CRSP Subsample with Daily Stock and Event Variables
+# Restrictions will be applied later
+# Select variables from the CRSP daily stock and event datasets
+crsp = conn.raw_sql("""
+                      select a.prc, a.ret, a.retx, a.shrout, a.vol, a.cfacpr, a.cfacshr, a.date, a.permno, a.permco,
+                      b.ticker, b.ncusip, b.shrcd, b.exchcd
+                      from crsp.dsf as a
+                      left join crsp.dsenames as b
+                      on a.permno=b.permno
+                      and b.namedt<=a.date
+                      and a.date<=b.nameendt
+                      where a.date >= '01/01/1959'
+                      and b.exchcd between 1 and 3
+                      """)
+
+# change variable format to int
+crsp[['permco', 'permno', 'shrcd', 'exchcd']] = crsp[['permco', 'permno', 'shrcd', 'exchcd']].astype(int)
+
+# Line up date to be end of month
+crsp['date'] = pd.to_datetime(crsp['date'])
+crsp['monthend'] = crsp['date'] + MonthEnd(0)  # set all the date to the standard end date of month
+
+crsp['me'] = crsp['prc'].abs() * crsp['shrout']  # calculate market equity
+
+# if Market Equity is Nan then let return equals to 0
+crsp['ret'] = np.where(crsp['me'].isnull(), 0, crsp['ret'])
+crsp['retx'] = np.where(crsp['me'].isnull(), 0, crsp['retx'])
+
+# impute me
+crsp = crsp.sort_values(by=['permno', 'date']).drop_duplicates()
+crsp['me'] = np.where(crsp['permno'] == crsp['permno'].shift(1), crsp['me'].fillna(method='ffill'), crsp['me'])
+
+# Aggregate Market Cap
+'''
+There are cases when the same firm (permco) has two or more securities (permno) at same date.
+For the purpose of ME for the firm, we aggregated all ME for a given permco, date.
+This aggregated ME will be assigned to the permno with the largest ME.
+'''
+# sum of me across different permno belonging to same permco a given date
+crsp_summe = crsp.groupby(['monthend', 'permco'])['me'].sum().reset_index()
+# largest mktcap within a permco/date
+crsp_maxme = crsp.groupby(['monthend', 'permco'])['me'].max().reset_index()
+# join by monthend/maxme to find the permno
+crsp1 = pd.merge(crsp, crsp_maxme, how='inner', on=['monthend', 'permco', 'me'])
+# drop me column and replace with the sum me
+crsp1 = crsp1.drop(['me'], axis=1)
+# join with sum of me to get the correct market cap info
+crsp2 = pd.merge(crsp1, crsp_summe, how='inner', on=['monthend', 'permco'])
+# sort by permno and date and also drop duplicates
+crsp2 = crsp2.sort_values(by=['permno', 'monthend']).drop_duplicates()
+
+
+#######################################################################################################################
+#                                                  Calculate                                                    #
+#######################################################################################################################
+
+
+def mom_1(start, end, df):
+    """
+    :param start: Order of starting lag
+    :param end: Order of ending lag
+    :param df: Dataframe
+    :return: Momentum factor
+    """
+    lag = pd.DataFrame()
+    result = 0
+    for i in range(start, end):
+        lag['mom%s' % i] = df.groupby(['permno'])['dtvm'].shift(i)
+        result = result + (lag['mom%s' % i])
+    result = result/(end-start)
+    return result
+
+
+def mom_2(start, end, df):
+    """
+    :param start: Order of starting lag
+    :param end: Order of ending lag
+    :param df: Dataframe
+    :return: Momentum factor
+    """
+    lag = pd.DataFrame()
+    result = 0
+    for i in range(start, end):
+        lag['mom%s' % i] = df.groupby(['permno'])['count'].shift(i)
+        result = result + (lag['mom%s' % i])
+    result = result
+    return result
+
+
+#Calculating daily trading volume
+#at least there is only one datapoint for a permno at one day
+crsp2['dtv'] = crsp2['prc'] * crsp2['vol']
+
+#Average trading value for one month
+#This will be a dataframe with every permno's average dtv for every month
+#It has 3 columns, and we reset its index
+dtv_m = crsp2.groupby(['permno','monthend'])[['dtv']].mean()
+dtv_m.rename(columns = {'dtv':'dtvm'})
+
+#record how many datapoints we have for a typical month
+dtv_m['month_count'] = crsp2.groupby(['permno','monthend'])[['dtv']].count()['dtv']
+dtv_m.reset_index()
+
+#Merge it back with crsp2 to only store monthly average data
+crsp3 = pd.merge(crsp2,dtv_m,how = 'inner', on = ['monthend','permno'])
+crsp3.drop(['dtv'],axis = 1)
+#Then calculate different dtvs
+
+#dtv
+crsp3['half_year_count'] =  mom_2(crsp3,0,6)
+crsp3['dtv'] = mom_1(crsp3,0,6)
+
+#change the ones with less than 50 records to nan
+crsp3['dtv'] = np.where(crsp3['half_year_count']<50, np.nan, crsp3['dtv'])
+
+with open('dtv.pkl', 'wb') as f:
+    pkl.dump(crsp3, f)
\ No newline at end of file

From 2473ad0a600159a42057d83361a276103a8eed9b Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Sat, 26 Sep 2020 17:47:22 +0800
Subject: [PATCH 24/28] Update dtv.py

---
 pychars/dtv.py | 102 +++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 94 insertions(+), 8 deletions(-)

diff --git a/pychars/dtv.py b/pychars/dtv.py
index 1f0e409..ca3274a 100644
--- a/pychars/dtv.py
+++ b/pychars/dtv.py
@@ -4,7 +4,9 @@
 import wrds
 from dateutil.relativedelta import *
 from pandas.tseries.offsets import *
+import datetime
 import pickle as pkl
+import multiprocessing as mp
 ###################
 # Connect to WRDS #
 ###################
@@ -81,7 +83,7 @@ def mom_1(start, end, df):
     lag = pd.DataFrame()
     result = 0
     for i in range(start, end):
-        lag['mom%s' % i] = df.groupby(['permno'])['dtvm'].shift(i)
+        lag['mom%s' % i] = df['dtvm'].shift(i)
         result = result + (lag['mom%s' % i])
     result = result/(end-start)
     return result
@@ -97,7 +99,7 @@ def mom_2(start, end, df):
     lag = pd.DataFrame()
     result = 0
     for i in range(start, end):
-        lag['mom%s' % i] = df.groupby(['permno'])['count'].shift(i)
+        lag['mom%s' % i] = df.groupby(['permno'])['day_count'].shift(i)
         result = result + (lag['mom%s' % i])
     result = result
     return result
@@ -114,20 +116,104 @@ def mom_2(start, end, df):
 dtv_m.rename(columns = {'dtv':'dtvm'})
 
 #record how many datapoints we have for a typical month
-dtv_m['month_count'] = crsp2.groupby(['permno','monthend'])[['dtv']].count()['dtv']
+dtv_m['day_count'] = crsp2.groupby(['permno','monthend'])[['dtv']].count()['dtv']
 dtv_m.reset_index()
 
 #Merge it back with crsp2 to only store monthly average data
 crsp3 = pd.merge(crsp2,dtv_m,how = 'inner', on = ['monthend','permno'])
 crsp3.drop(['dtv'],axis = 1)
-#Then calculate different dtvs
 
-#dtv
-crsp3['half_year_count'] =  mom_2(crsp3,0,6)
-crsp3['dtv'] = mom_1(crsp3,0,6)
+#Generate firm list
+df_firm = crsp3.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+#Extract number of data points for each permno
+crsp3['month_count'] = crsp3.groupby('permno').cumcount()
+month_num = crsp3.groupby('permno')['month_count'].tail(1)
+month_num = month_num.astype(int)
 
+#dtv
+#crsp3['half_year_count'] =  mom_2(crsp3,0,6)
+#crsp3['dtv'] = mom_1(crsp3,0,6)
 #change the ones with less than 50 records to nan
-crsp3['dtv'] = np.where(crsp3['half_year_count']<50, np.nan, crsp3['dtv'])
+#crsp3['dtv'] = np.where(crsp3['half_year_count']<50, np.nan, crsp3['dtv'])
+
+
+def get_dtv(df, firm_list):
+    """
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            temp = df[(df['permno'] == firm) & (i - 5 <= df['month_count']) & (df['month_count'] <= i)]
+            # if observations in last 3 months are less than 2 months, we drop the rvar of this month
+            if temp['permno'].count() < 2:
+                pass
+            else:
+                index = temp.tail(1).index
+                df.loc[index, 'dtv'] = mom_1(0,6,temp)
+    return df
+
+
+def sub_df(start, end, step):
+    """
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['cr' + str(h)] = pd.merge(crsp3, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['cr' + str(h)] = pd.merge(crsp3, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+def main(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_dtv, (df['cr%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+if __name__ == '__main__':
+    crsp3 = main(0, 1, 0.05)
+
+crsp3 = crsp3.dropna(subset=['dtv'])  # drop NA due to rolling
+crsp3 = crsp3.reset_index(drop=True)
+crsp3 = crsp3[['permno', 'date', 'dtv']]
+
 
 with open('dtv.pkl', 'wb') as f:
     pkl.dump(crsp3, f)
\ No newline at end of file

From 8beac9b4dcc4283c2837fe9cb2135d7c97d97f48 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Sun, 27 Sep 2020 18:19:48 +0800
Subject: [PATCH 25/28] Update friction chars

---
 pychars/hxz_Isff.py |   0
 pychars/hxz_Ivff.py | 199 ++++++++++++++++++++++++++++++++++++++++++++
 pychars/hxz_Ivq.py  |   0
 pychars/hxz_tv.py   | 153 ++++++++++++++++++++++++++++++++++
 4 files changed, 352 insertions(+)
 create mode 100644 pychars/hxz_Isff.py
 create mode 100644 pychars/hxz_Ivff.py
 create mode 100644 pychars/hxz_Ivq.py
 create mode 100644 pychars/hxz_tv.py

diff --git a/pychars/hxz_Isff.py b/pychars/hxz_Isff.py
new file mode 100644
index 0000000..e69de29
diff --git a/pychars/hxz_Ivff.py b/pychars/hxz_Ivff.py
new file mode 100644
index 0000000..0d37111
--- /dev/null
+++ b/pychars/hxz_Ivff.py
@@ -0,0 +1,199 @@
+# Fama & French 3 factors residual variance
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block
+crsp = conn.raw_sql("""
+                      select a.permno, a.date, a.ret, (a.ret - b.rf) as exret, b.mktrf, b.smb, b.hml
+                      from crsp.dsf as a
+                      left join ff.factors_daily as b
+                      on a.date=b.date
+                      where a.date > '01/01/1959'
+                      """)
+
+# sort variables by permno and date
+crsp = crsp.sort_values(by=['permno', 'date'])
+
+# change variable format to int
+crsp['permno'] = crsp['permno'].astype(int)
+
+# Line up date to be end of month
+crsp['date'] = pd.to_datetime(crsp['date'])
+
+# find the closest trading day to the end of the month
+crsp['monthend'] = crsp['date'] + MonthEnd(0)
+crsp['date_diff'] = crsp['monthend'] - crsp['date']
+date_temp = crsp.groupby(['permno', 'monthend'])['date_diff'].min()
+date_temp = pd.DataFrame(date_temp)  # convert Series to DataFrame
+date_temp.reset_index(inplace=True)
+date_temp.rename(columns={'date_diff': 'min_diff'}, inplace=True)
+crsp = pd.merge(crsp, date_temp, how='left', on=['permno', 'monthend'])
+crsp['sig'] = np.where(crsp['date_diff'] == crsp['min_diff'], 1, np.nan)
+
+# label every date of month end
+crsp['month_count'] = crsp[crsp['sig'] == 1].groupby(['permno']).cumcount()
+
+# label numbers of months for a firm
+month_num = crsp[crsp['sig'] == 1].groupby(['permno'])['month_count'].tail(1)
+month_num = month_num.astype(int)
+
+# mark the number of each month to each day of this month
+crsp['month_count'] = crsp.groupby(['permno'])['month_count'].fillna(method='bfill')
+
+# crate a firm list
+df_firm = crsp.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+######################
+# Calculate the beta #
+######################
+# function that get multiple beta
+''''
+rolling_window = 60  # 60 trading days
+crsp['beta_mktrf'] = np.nan
+crsp['beta_smb'] = np.nan
+crsp['beta_hml'] = np.nan
+
+
+def get_beta(df):
+    """
+    The original idea of calculate beta is using formula (X'MX)^(-1)X'MY,
+    where M = I - 1(1'1)^{-1}1, I is a identity matrix.
+
+    """
+    temp = crsp.loc[df.index]  # extract the rolling sub dataframe from original dataframe
+    X = np.mat(temp[['mktrf', 'smb', 'hml']])
+    Y = np.mat(temp[['exret']])
+    ones = np.mat(np.ones(rolling_window)).T
+    M = np.identity(rolling_window) - ones.dot((ones.T.dot(ones)).I).dot(ones.T)
+    beta = (X.T.dot(M).dot(X)).I.dot((X.T.dot(M).dot(Y)))
+    crsp['beta_mktrf'].loc[df.index[-1:]] = beta[0]
+    crsp['beta_smb'].loc[df.index[-1:]] = beta[1]
+    crsp['beta_hml'].loc[df.index[-1:]] = beta[2]
+    return 0  # we do not need the rolling outcome since rolling cannot return different values in different columns
+
+
+# calculate beta through rolling window
+crsp_temp = crsp.groupby('permno').rolling(rolling_window).apply(get_beta, raw=False)
+'''
+
+######################
+# Calculate residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    """
+
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            temp = df[(df['permno'] == firm) & (i - 1 <= df['month_count']) & (df['month_count'] < i)]
+            # if observations in last 3 months are less 21, we drop the rvar of this month
+            if temp['permno'].count() < 21:
+                pass
+            else:
+                rolling_window = temp['permno'].count()
+                index = temp.tail(1).index
+                X = pd.DataFrame()
+                X[['mktrf', 'smb', 'hml']] = temp[['mktrf', 'smb', 'hml']]
+                X['intercept'] = 1
+                X = X[['intercept', 'mktrf', 'smb', 'hml']]
+                X = np.mat(X)
+                Y = np.mat(temp[['exret']])
+                res = (np.identity(rolling_window) - X.dot(X.T.dot(X).I).dot(X.T)).dot(Y)
+                res_var = res.var(ddof=1)
+                df.loc[index, 'rvar'] = res_var
+    return df
+
+
+def sub_df(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+
+def main(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_res_var, (df['crsp%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+# process dataframe
+crsp = crsp.dropna(subset=['rvar'])  # drop NA due to rolling
+crsp = crsp.rename(columns={'rvar': 'Ivff'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'Ivff']]
+
+with open('Ivff.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
\ No newline at end of file
diff --git a/pychars/hxz_Ivq.py b/pychars/hxz_Ivq.py
new file mode 100644
index 0000000..e69de29
diff --git a/pychars/hxz_tv.py b/pychars/hxz_tv.py
new file mode 100644
index 0000000..a15d04d
--- /dev/null
+++ b/pychars/hxz_tv.py
@@ -0,0 +1,153 @@
+# CAPM residual variance
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block
+crsp = conn.raw_sql("""
+                      select permno, date, ret
+                      from crsp.dsf
+                      """)
+
+# sort variables by permno and date
+crsp = crsp.sort_values(by=['permno', 'date'])
+
+# change variable format to int
+crsp['permno'] = crsp['permno'].astype(int)
+
+# Line up date to be end of month
+crsp['date'] = pd.to_datetime(crsp['date'])
+
+# find the closest trading day to the end of the month
+crsp['monthend'] = crsp['date'] + MonthEnd(0)
+crsp['date_diff'] = crsp['monthend'] - crsp['date']
+date_temp = crsp.groupby(['permno', 'monthend'])['date_diff'].min()
+date_temp = pd.DataFrame(date_temp)  # convert Series to DataFrame
+date_temp.reset_index(inplace=True)
+date_temp.rename(columns={'date_diff': 'min_diff'}, inplace=True)
+crsp = pd.merge(crsp, date_temp, how='left', on=['permno', 'monthend'])
+crsp['sig'] = np.where(crsp['date_diff'] == crsp['min_diff'], 1, np.nan)
+
+# label every date of month end
+crsp['month_count'] = crsp[crsp['sig'] == 1].groupby(['permno']).cumcount()
+# label numbers of months for a firm
+month_num = crsp[crsp['sig'] == 1].groupby(['permno'])['month_count'].tail(1)
+month_num = month_num.astype(int)
+
+# mark the number of each month to each day of this month
+crsp['month_count'] = crsp.groupby(['permno'])['month_count'].fillna(method='bfill')
+
+# crate a firm list
+df_firm = crsp.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+######################
+# Calculate residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    """
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            temp = df[(df['permno'] == firm) & (i - 1 <= df['month_count']) & (df['month_count'] < i)]
+            # if observations in last 3 months are less 21, we drop the rvar of this month
+            if temp['permno'].count() < 15:
+                pass
+            else:
+                index = temp.tail(1).index
+                rvar = temp['ret'].var(ddof = 1)
+                df.loc[index, 'rvar'] = rvar
+    return df
+
+
+def sub_df(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+
+def main(start, end, step):
+    """
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_res_var, (df['crsp%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+# process dataframe
+crsp = crsp.dropna(subset=['rvar'])  # drop NA due to rolling
+crsp = crsp.rename(columns={'rvar': 'tv'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'rvar_capm']]
+
+with open('tv.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
\ No newline at end of file

From be583a41aebce9177227b82a3c277d432ca36170 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Mon, 28 Sep 2020 00:32:10 +0800
Subject: [PATCH 26/28] Update hxz_Ivq.py

---
 pychars/hxz_Ivq.py | 176 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 176 insertions(+)

diff --git a/pychars/hxz_Ivq.py b/pychars/hxz_Ivq.py
index e69de29..c659263 100644
--- a/pychars/hxz_Ivq.py
+++ b/pychars/hxz_Ivq.py
@@ -0,0 +1,176 @@
+# CAPM residual variance
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block
+crsp = conn.raw_sql("""
+                      select permno, date, ret
+                      from crsp.dsf
+                      where date >= '01/03/1967'
+                      """)
+
+
+qmodel = pd.read_csv("q5_factors_daily_2019a.csv")
+
+# sort variables by permno and date
+crsp = crsp.sort_values(by=['permno', 'date'])
+
+# change variable format to int
+crsp['permno'] = crsp['permno'].astype(int)
+
+# Line up date to be end of month
+# Merge crsp and q factor
+crsp['date'] = pd.to_datetime(crsp['date'])
+qmodel.rename(columns = {'DATE':'date'}, inplace = True)
+qmodel['date'] = pd.to_datetime(qmodel['date'],format='%Y%m%d', errors='ignore')
+print(qmodel)
+crsp = pd.merge(crsp,qmodel,how = 'inner', on = ['date'])
+print(crsp)
+
+# find the closest trading day to the end of the month
+crsp['monthend'] = crsp['date'] + MonthEnd(0)
+crsp['date_diff'] = crsp['monthend'] - crsp['date']
+date_temp = crsp.groupby(['permno', 'monthend'])['date_diff'].min()
+date_temp = pd.DataFrame(date_temp)  # convert Series to DataFrame
+date_temp.reset_index(inplace=True)
+date_temp.rename(columns={'date_diff': 'min_diff'}, inplace=True)
+crsp = pd.merge(crsp, date_temp, how='left', on=['permno', 'monthend'])
+crsp['sig'] = np.where(crsp['date_diff'] == crsp['min_diff'], 1, np.nan)
+
+# label every date of month end
+crsp['month_count'] = crsp[crsp['sig'] == 1].groupby(['permno']).cumcount()
+# label numbers of months for a firm
+month_num = crsp[crsp['sig'] == 1].groupby(['permno'])['month_count'].tail(1)
+month_num = month_num.astype(int)
+
+# mark the number of each month to each day of this month
+crsp['month_count'] = crsp.groupby(['permno'])['month_count'].fillna(method='bfill')
+
+# crate a firm list
+df_firm = crsp.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+######################
+# Calculate residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    """
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            temp = df[(df['permno'] == firm) & (i - 1 <= df['month_count']) & (df['month_count'] < i)]
+            # if observations in last 1 month are less 15, we drop the rvar of this month
+            if temp['permno'].count() < 15:
+                pass
+            else:
+                rolling_window = temp['permno'].count()
+                index = temp.tail(1).index
+                X = pd.DataFrame()
+                temp['rmf'] = temp['R_MKT'] - temp['R_F']
+                temp['y'] = temp['ret'] - temp['R_F']
+                X[['rme']] = temp[['R_ME']]
+                X[['rm_f']] = temp[['rmf']]
+                X[['ria']] = temp[['R_IA']]
+                X[['roe']] = temp[['R_ROE']]
+                X['intercept'] = 1
+                X = X[['intercept', 'rm_f','rme','ria','roe']]
+                X = np.mat(X)
+                Y = np.mat(temp[['y']])
+                res = (np.identity(rolling_window) - X.dot(X.T.dot(X).I).dot(X.T)).dot(Y)
+                res_var = res.var(ddof=1)
+                df.loc[index, 'rvar'] = res_var
+    return df
+
+
+def sub_df(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+
+def main(start, end, step):
+    """
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_res_var, (df['crsp%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+# process dataframe
+crsp = crsp.dropna(subset=['rvar'])  # drop NA due to rolling
+crsp = crsp.rename(columns={'rvar': 'Ivq'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'Ivq']]
+
+with open('Ivq.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
\ No newline at end of file

From 9ed0389f6639830bdaeeb499a74503344927810a Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Mon, 28 Sep 2020 11:04:13 +0800
Subject: [PATCH 27/28] Create hxz_sv.py

---
 pychars/hxz_sv.py | 174 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 174 insertions(+)
 create mode 100644 pychars/hxz_sv.py

diff --git a/pychars/hxz_sv.py b/pychars/hxz_sv.py
new file mode 100644
index 0000000..7a725d4
--- /dev/null
+++ b/pychars/hxz_sv.py
@@ -0,0 +1,174 @@
+# CAPM residual variance
+# Note: Please use the latest version of pandas, this version should support returning to pd.Series after rolling
+# To get a faster speed, we split the big dataframe into small ones
+# Then using different process to calculate the variance
+# We use 20 process to calculate variance, you can change the number of process according to your CPU situation
+# You can use the following code to check your CPU situation
+# import multiprocessing
+# multiprocessing.cpu_count()
+
+import pandas as pd
+import numpy as np
+import datetime as dt
+import wrds
+from dateutil.relativedelta import *
+from pandas.tseries.offsets import *
+import datetime
+import pickle as pkl
+import multiprocessing as mp
+
+###################
+# Connect to WRDS #
+###################
+conn = wrds.Connection()
+
+# CRSP Block
+crsp = conn.raw_sql("""
+                      select a.permno, a.date, a.ret, (a.ret - b.rf) as exret, b.mktrf
+                      from crsp.dsf as a
+                      left join ff.factors_daily as b
+                      on a.date=b.date
+                      where a.date >= '01/01/1986'
+                      """)
+
+cboe = conn.raw_sql("""
+                      select vxoh, vxol
+                      from cboe
+                      where date >= '01/01/1986'
+                      """)
+
+# sort variables by permno and date
+crsp = crsp.sort_values(by=['permno', 'date'])
+
+# change variable format to int
+crsp['permno'] = crsp['permno'].astype(int)
+
+# Line up date to be end of month
+crsp['date'] = pd.to_datetime(crsp['date'])
+cboe['date'] = pd.to_datetime(cboe['date'])
+crsp = pd.merge(crsp,cboe,how = 'left',on = ['date'])
+
+# find the closest trading day to the end of the month
+crsp['monthend'] = crsp['date'] + MonthEnd(0)
+crsp['date_diff'] = crsp['monthend'] - crsp['date']
+date_temp = crsp.groupby(['permno', 'monthend'])['date_diff'].min()
+date_temp = pd.DataFrame(date_temp)  # convert Series to DataFrame
+date_temp.reset_index(inplace=True)
+date_temp.rename(columns={'date_diff': 'min_diff'}, inplace=True)
+crsp = pd.merge(crsp, date_temp, how='left', on=['permno', 'monthend'])
+crsp['sig'] = np.where(crsp['date_diff'] == crsp['min_diff'], 1, np.nan)
+
+# label every date of month end
+crsp['month_count'] = crsp[crsp['sig'] == 1].groupby(['permno']).cumcount()
+# label numbers of months for a firm
+month_num = crsp[crsp['sig'] == 1].groupby(['permno'])['month_count'].tail(1)
+month_num = month_num.astype(int)
+
+# mark the number of each month to each day of this month
+crsp['month_count'] = crsp.groupby(['permno'])['month_count'].fillna(method='bfill')
+
+# crate a firm list
+df_firm = crsp.drop_duplicates(['permno'])
+df_firm = df_firm[['permno']]
+df_firm['permno'] = df_firm['permno'].astype(int)
+df_firm = df_firm.reset_index(drop=True)
+df_firm = df_firm.reset_index()
+df_firm = df_firm.rename(columns={'index': 'count'})
+
+######################
+# Calculate residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    """
+
+    :param df: stock dataframe
+    :param firm_list: list of firms matching stock dataframe
+    :return: dataframe with variance of residual
+    """
+    for firm, count, prog in zip(firm_list['permno'], month_num, range(firm_list['permno'].count()+1)):
+        prog = prog + 1
+        print('processing permno %s' % firm, '/', 'finished', '%.2f%%' % ((prog/firm_list['permno'].count())*100))
+        for i in range(count + 1):
+            # if you want to change the rolling window, please change here: i - 2 means 3 months is a window.
+            temp = df[(df['permno'] == firm) & (i - 1 <= df['month_count']) & (df['month_count'] < i)]
+            # if observations in last 3 months are less 21, we drop the rvar of this month
+            if temp['permno'].count() < 15:
+                pass
+            else:
+                temp['dvxo'] = temp['vxoh'] - temp['vxol']
+                rolling_window = temp['permno'].count()
+                index = temp.tail(1).index
+                X = pd.DataFrame()
+                X[['mktrf']] = temp[['mktrf']]
+                X[['dvxo']] = temp['dvxo']
+                X['intercept'] = 1
+                X = X[['intercept', 'mktrf', 'dvxo']]
+                X = np.mat(X)
+                Y = np.mat(temp[['exret']])
+                bt = (X.T.dot(X).I).dot(X.T).dot(Y)
+                df.loc[index, 'rvar'] = np.float(bt[2])
+    return df
+
+
+def sub_df(start, end, step):
+    """
+
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dictionary including all the 'firm_list' dataframe and 'stock data' dataframe
+    """
+    # we use dict to store different sub dataframe
+    temp = {}
+    for i, h in zip(np.arange(start, end, step), range(int((end-start)/step))):
+        print('processing splitting dataframe:', round(i, 2), 'to', round(i + step, 2))
+        if i == 0:  # to get the left point
+            temp['firm' + str(h)] = df_firm[df_firm['count'] <= df_firm['count'].quantile(i + step)]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+        else:
+            temp['firm' + str(h)] = df_firm[(df_firm['count'].quantile(i) < df_firm['count']) & (
+                    df_firm['count'] <= df_firm['count'].quantile(i + step))]
+            temp['crsp' + str(h)] = pd.merge(crsp, temp['firm' + str(h)], how='left',
+                                             on='permno').dropna(subset=['count'])
+    return temp
+
+
+def main(start, end, step):
+    """
+    :param start: the quantile to start cutting, usually it should be 0
+    :param end: the quantile to end cutting, usually it should be 1
+    :param step: quantile step
+    :return: a dataframe with calculated variance of residual
+    """
+    df = sub_df(start, end, step)
+    pool = mp.Pool()
+    p_dict = {}
+    for i in range(int((end-start)/step)):
+        p_dict['p' + str(i)] = pool.apply_async(get_res_var, (df['crsp%s' % i], df['firm%s' % i],))
+    pool.close()
+    pool.join()
+    result = pd.DataFrame()
+    print('processing pd.concat')
+    for h in range(int((end-start)/step)):
+        result = pd.concat([result, p_dict['p%s' % h].get()])
+    return result
+
+
+
+# calculate variance of residual through rolling window
+# Note: please split dataframe according to your CPU situation. For example, we split dataframe to (1-0)/0.05 = 20 sub
+# dataframes here, so the function will use 20 cores to calculate variance of residual.
+if __name__ == '__main__':
+    crsp = main(0, 1, 0.05)
+
+# process dataframe
+crsp = crsp.dropna(subset=['rvar'])  # drop NA due to rolling
+crsp = crsp.rename(columns={'rvar': 'sv'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'sv']]
+
+with open('sv.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
\ No newline at end of file

From 780873c7698f0ced708dfff4d655dfc9ce07ea41 Mon Sep 17 00:00:00 2001
From: TianXie1999 <55803615+TianXie1999@users.noreply.github.com>
Date: Fri, 9 Oct 2020 21:37:18 +0800
Subject: [PATCH 28/28] Update hxz_Ivq.py

---
 pychars/hxz_Ivq.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pychars/hxz_Ivq.py b/pychars/hxz_Ivq.py
index c659263..ec21d51 100644
--- a/pychars/hxz_Ivq.py
+++ b/pychars/hxz_Ivq.py
@@ -29,7 +29,7 @@
                       where date >= '01/03/1967'
                       """)
 
-
+#Download data from http://global-q.org/factors.html
 qmodel = pd.read_csv("q5_factors_daily_2019a.csv")
 
 # sort variables by permno and date
@@ -43,9 +43,8 @@
 crsp['date'] = pd.to_datetime(crsp['date'])
 qmodel.rename(columns = {'DATE':'date'}, inplace = True)
 qmodel['date'] = pd.to_datetime(qmodel['date'],format='%Y%m%d', errors='ignore')
-print(qmodel)
 crsp = pd.merge(crsp,qmodel,how = 'inner', on = ['date'])
-print(crsp)
+
 
 # find the closest trading day to the end of the month
 crsp['monthend'] = crsp['date'] + MonthEnd(0)