diff --git a/pychars/accounting.py b/pychars/accounting.py
index 41e7194..d252b85 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 #
 ###################
@@ -59,20 +58,20 @@ 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,
                     
                     /*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,
+                    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_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
@@ -191,7 +190,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)
 
@@ -199,6 +198,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'])
 
@@ -210,8 +210,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
 
 # count single stock years
 data_rawa['count'] = data_rawa.groupby(['gvkey']).cumcount()
@@ -270,10 +270,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，checked
 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)
+
+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)
@@ -282,7 +320,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'])
@@ -292,6 +330,16 @@ 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, 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, np.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']
@@ -306,7 +354,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']
 
@@ -316,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
@@ -371,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'])
@@ -433,6 +481,75 @@ 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['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)
+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'])
+#data_rawa['dwc'] = data_rawa['dwc']/data_rawa['at_l1']
+
+#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['2ig'] = 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'] 
+
+#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'])
+
+#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']
+
+#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',
                      'sic', 'acc', 'agr', 'bm', 'cfp', 'ep', 'ni', 'op', 'rsup', 'cash', 'chcsho',
@@ -502,8 +619,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
 
 # count single stock years
 data_rawq['count'] = data_rawq.groupby(['gvkey']).cumcount()
@@ -572,7 +689,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']
 
@@ -630,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'])
@@ -714,6 +834,21 @@ 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
+
+#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']
+
+#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',
@@ -734,6 +869,77 @@ 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
+
+
+
+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 
@@ -817,4 +1023,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)
+
+
diff --git a/pychars/dtv.py b/pychars/dtv.py
new file mode 100644
index 0000000..ca3274a
--- /dev/null
+++ b/pychars/dtv.py
@@ -0,0 +1,219 @@
+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                                                    #
+#######################################################################################################################
+
+# 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['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'])['day_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['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)
+
+#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'])
+
+
+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
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..ec21d51
--- /dev/null
+++ b/pychars/hxz_Ivq.py
@@ -0,0 +1,175 @@
+# 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'
+                      """)
+
+#Download data from http://global-q.org/factors.html
+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')
+crsp = pd.merge(crsp,qmodel,how = 'inner', 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 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
diff --git a/pychars/hxz_dtv.py b/pychars/hxz_dtv.py
new file mode 100644
index 0000000..f678b2e
--- /dev/null
+++ b/pychars/hxz_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
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
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
diff --git a/pychars/rmom_capm.py b/pychars/rmom_capm.py
new file mode 100644
index 0000000..fcdc706
--- /dev/null
+++ b/pychars/rmom_capm.py
@@ -0,0 +1,155 @@
+# 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
+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 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']]
+
+# 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 residual #
+######################
+
+
+def get_res_var(df, firm_list):
+    #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() < 20:
+                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)
+                #print(df)
+                df.loc[index, 'res'] = res[-1]
+    return df
+
+
+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))):
+        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):
+    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_capm_1m'})
+crsp = crsp.reset_index(drop=True)
+crsp = crsp[['permno', 'date', 'rmom_capm_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_capm_1m'].shift(i)
+        result = result * (1+lag['mom%s' % i])
+    result = result - 1
+    return result
+
+
+crsp['rmom_capm_12m'] = mom(1,12,crsp)
+crsp['rmom_capm_60m'] = mom(12,60,crsp)
+
+with open('rmom_capm.pkl', 'wb') as f:
+    pkl.dump(crsp, f)
+
+
+
+
+
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
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