<jupyter_start><jupyter_text>Part 1: For each indices we calibrate the model's parameters an draw mean levels against market caps<jupyter_code># read the data to a Dataframe
df_cap = pd.read_csv('dja_cap.csv')
# Create the Scaler object
#scaler = preprocessing.StandardScaler()
# Fit your data on the scaler object
#df_cap = pd.DataFrame(scaler.fit_transform(df_cap))
# add dates
dates = pd.bdate_range(start='2010-01-04', end=None, periods=df_cap.shape[0], freq='B')
df_cap['date'] = dates
df_cap.set_index('date',inplace=True)

# get stock indices
stock_id=df_cap.iloc[0,:].index



df_cap.head()<jupyter_output><empty_output><jupyter_text>Perform Parameter's calibrationAssume $q_t = log X_t$ and apply Ito's lemma:$$\Delta q_t = \kappa_x \cdot \left( Wz^{'}_{t} - \frac{1}{2} \frac{\sigma^2}{\kappa_x} - x_t\right) \Delta t + \sigma \cdot \Delta W_t$$The problem reduces to find that set of parameters so that, $$  \underset{\Theta}{Arg Min}  \; \sum_{t=0}^{T-1} \left| \Sigma_x \right|  +  \left( {\bf v}_t \right)^{T}  \Sigma_x^{-1}   {\bf v}_t + \lambda (w1 + w2 - 1)^2$$where$$  {\bf v}_t \equiv log \frac{{\bf x}_{t+1} }{{\bf x}_{t}}  -  \kappa_x \circ \left({\bf W} {\bf z}_t' - \frac{1}{2} \frac{\sigma^2}{\kappa_x}  - {\bf x}_t \right)  $$<jupyter_code># likelihood function definition
def mle(params,X,Z1,Z2,asset_id):
    kappa = params[0]
    sigma = params[1]
    correl = 0.0

    w0 = params[2]
    w1 = params[3]
    w2 = params[4]
    X_t = np.array(X.iloc[0:-1,:])    
    Z1 = np.array(Z1.iloc[:-1, :])
    Z2 = np.array(Z2.iloc[:-1, :])

   
    log_returns_X = np.array(np.log(1 + X.pct_change()).iloc[1:, :]) # removed first NaN




    nb_assets = log_returns_X.shape[1]
    sample_size = log_returns_X.shape[0]

        
    # WZ
    Wz = w0 + w1 * Z1 + w2 * Z2
 
    # VCV Matrix inverse ...
    Sigma = np.eye(log_returns_X.shape[1]) * (sigma ** (2))


        # Calculate v_t
    v_t = log_returns_X - kappa * (Wz - 0.5 * sigma**2/kappa - X_t)
###################
    v_t=v_t[:,asset_id]
    nb_assets=1
####################
#print(np.dot(v_t[0], v_t[0]))    
        # Calculate v'v/sigma^2 (Valid as long as we define VCV to be diagonal => independence)                
    quadratic_term = 0
    for i in range(sample_size):            
            #density *= multivariate_normal.pdf(self.log_returns_X[i], mu_t[i], Sigma)
        quadratic_term += np.dot(v_t[i], v_t[i])
    quadratic_term /= (sigma**2)
   


        # Calculate objective function Sum |Sigma| + v'v/sigma^2
    objective_funct_val = quadratic_term + sample_size * ( nb_assets) * np.log(sigma)

        # Calculate penalty function
    penalty_coef =  10e-2;
    penalty_function =  penalty_coef * (w1 + w2 - 1)**2
    return objective_funct_val + penalty_function
def print_estimated_mean_rate(date_from, date_to, asset_id, parameters):

    # Retrieve market signals ...
    Z1_array = np.array(Z1.loc[date_from:date_to, asset_id])
    Z2_array = np.array(Z2.loc[date_from:date_to, asset_id])

    # Retrieve dates and mkt and model caps ..
    dates = np.array(X.loc[date_from:date_to, :].index)
    X_t_model = np.array(X.loc[date_from:date_to, asset_id])
    X_t_mket = np.array(X.loc[date_from:date_to, asset_id])

    
    kappa = parameters[0]
    sigma = parameters[1]
    w0 = parameters[2]
    w1 = parameters[3]
    w2 = parameters[4]
    
    # -> Calculate Wz
    Wz = w0 + w1 * Z1_array + w2 * Z2_array

    # lon run mean level 
    theta_z = Wz - 0.5 * sigma**2 / kappa
 
    


    # GRAPH ...
    fig = plt.figure(figsize=(15,6))

    # Print theta z
    ax = fig.add_subplot(1,2,1)
    ax.plot( theta_z, label = 'Theta')
    ax.legend(loc='best')
    ax.set_title('Mean rate')

    # print market caps ...
    ax = fig.add_subplot(1,2,2)
    ax.plot( X_t_mket, label = 'Mkt cap ')
    ax.plot( theta_z, label = 'Mean rate')
    ax.legend(loc='best')
    ax.set_title('Normalized capiltalization')


    plt.show()<jupyter_output><empty_output><jupyter_text>Example of calibration and check for 'AXP'  stock<jupyter_code>date_from = '2011-01-01'
date_to = '2012-09-01' 


# start data analisis and signal extraction

X = df_cap.loc[:,:]
X_mean = X.mean(axis = 0).sum() 
X /= X_mean 

window_1 = 15  
window_2 = 30
         # Equal weighted
Z1 = X.rolling(window=window_1).mean()
Z2 = X.rolling(window=window_2).mean()

# remove NaN
X = X.iloc[window_2:,:]
Z1 = Z1.iloc[window_2:,:]
Z2 = Z2.iloc[window_2:,:]
X_window = X.loc[date_from:date_to, :]
Z1_window = Z1.loc[date_from:date_to, :]
Z2_window = Z2.loc[date_from:date_to, :]

params=[0.2,0.2,0.2,0.2,0.2]
lb = np.array([-10, 0.0001, -5, 0.0001, 0.0001])
ub = np.array([10,   5,     5, 1,      1])
bounds = [(lb[0], ub[0]), (lb[1], ub[1]), (lb[2], ub[2]), (lb[3], ub[3]), (lb[4], ub[4])]
# PERFORM CALIBRATION
asset_id = 1 # 1 corresponds to 'AXP'
res = minimize(mle,params,args=(X_window,Z1_window,Z2_window,asset_id,),bounds = bounds,
                       method='L-BFGS-B', options={'gtol': 1e-10, 'disp': True}, tol = 1e-9)
print(res)

# check 
print_estimated_mean_rate(date_from, date_to, 'AXP', res.x)<jupyter_output>fun: -1406.4563541479683
 hess_inv: <5x5 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 0.0347427 , -0.04615686,  1.60034688,  0.62620984, -0.80319751])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 570
      nit: 62
   status: 0
  success: True
        x: array([4.0156735 , 0.02373729, 0.0042492 , 0.25963029, 0.41657371])<jupyter_text>calibration of the parameters for all the stock<jupyter_code>dfObj = pd.DataFrame(columns=['k', 'sigma', 'W0','W1','W2'],index=stock_id)
for asset_id in range(0,30):
    res = minimize(mle,params,args=(X_window,Z1_window,Z2_window,asset_id,),bounds = bounds,
                       method='L-BFGS-B', options={'gtol': 1e-10, 'disp': True}, tol = 1e-9)

# Append rows in Empty Dataframe by adding dictionaries
    dfObj.loc[stock_id[asset_id]]=[res.x[0],res.x[1],res.x[2],res.x[3],res.x[4]]
 
#print("Dataframe Contens ", dfObj, sep='\n')
dfObj<jupyter_output><empty_output><jupyter_text>Part 2: same as part I but with EWMA as a signal<jupyter_code>Z1 = X.ewm(alpha = 0.8, adjust = True).mean()
Z2 = X.ewm(alpha = 0.96, adjust = True).mean()
# remove NaN
X = X.iloc[window_2:,:]
Z1 = Z1.iloc[window_2:,:]
Z2 = Z2.iloc[window_2:,:]
date_from = '2011-01-01'
date_to = '2012-09-01' 

#-> Retrieve X, Z1 and Z2 only for the estimation window ...
X_window = X.loc[date_from:date_to, :]
Z1_window = Z1.loc[date_from:date_to, :]
Z2_window = Z2.loc[date_from:date_to, :]

params=[0.2,0.2,0.2,0.2,0.2]
lb = np.array([-0.01, 0.001, -5, 0.0001, 0.0001])
ub = np.array([5,   0.5,     5, 1,      1])
bounds = [(lb[0], ub[0]), (lb[1], ub[1]), (lb[2], ub[2]), (lb[3], ub[3]), (lb[4], ub[4])]
# PERFORM CALIBRATION
asset_id=1
res = minimize(mle,params,args=(X_window,Z1_window,Z2_window,asset_id,),bounds = bounds,
                       method='L-BFGS-B', options={'gtol': 1e-10, 'disp': True}, tol = 1e-9)
print(res)

print_estimated_mean_rate(date_from, date_to, 'AXP', res.x)
dfObj = pd.DataFrame(columns=['k', 'sigma', 'W0','W1','W2'],index=stock_id)
for asset_id in range(0,30):
    res = minimize(mle,params,args=(X_window,Z1_window,Z2_window,asset_id,),bounds = bounds,
                       method='L-BFGS-B', options={'gtol': 1e-10, 'disp': True}, tol = 1e-9)

# Append rows in Empty Dataframe by adding dictionaries
    dfObj.loc[stock_id[asset_id]]=[res.x[0],res.x[1],res.x[2],res.x[3],res.x[4]]
 
#print("Dataframe Contens ", dfObj, sep='\n')
dfObj<jupyter_output><empty_output><jupyter_text>Part III calibration of the parameters for S&P500:now we minimize the likelihood for all the stocks together<jupyter_code># likelihood function definition
def mleMarket(params,X,Z1,Z2):
    kappa = params[0]
    sigma = params[1]
    correl = 0.0

    w0 = params[2]
    w1 = params[3]
    w2 = params[4]
    X_t = np.array(X.iloc[0:-1,:])    
    Z1 = np.array(Z1.iloc[:-1, :])
    Z2 = np.array(Z2.iloc[:-1, :])

   
    log_returns_X = np.array(np.log(1 + X.pct_change()).iloc[1:, :]) # removed first NaN




    nb_assets = log_returns_X.shape[1]
    sample_size = log_returns_X.shape[0]

        
    # WZ
    Wz = w0 + w1 * Z1 + w2 * Z2
 
    # VCV Matrix inverse ...
    Sigma = np.eye(log_returns_X.shape[1]) * (sigma ** (2))


        # Calculate v_t
    v_t = log_returns_X - kappa * (Wz - 0.5 * sigma**2/kappa - X_t)
###################

####################
#print(np.dot(v_t[0], v_t[0]))    
        # Calculate v'v/sigma^2 (Valid as long as we define VCV to be diagonal => independence)                
    quadratic_term = 0
    for i in range(sample_size):            
            #density *= multivariate_normal.pdf(self.log_returns_X[i], mu_t[i], Sigma)
        quadratic_term += np.dot(v_t[i], v_t[i])
    quadratic_term /= (sigma**2)
   


        # Calculate objective function Sum |Sigma| + v'v/sigma^2
    objective_funct_val = quadratic_term + sample_size * ( nb_assets) * np.log(sigma)

        # Calculate penalty function
    penalty_coef =  10e-2;
    penalty_function =  penalty_coef * (w1 + w2 - 1)**2
    return objective_funct_val + penalty_function
# Put the # EUROSTOXX
df_cap = pd.read_csv('spx_holdings_and_spx_closeprice.csv')

dates = df_cap.iloc[:,0]

df_cap['date'] = dates

df_cap.set_index('date',inplace=True)
df_cap = df_cap.drop(df_cap.columns[[0]], axis=1)  # df.columns is zero-based pd.Index 
df_cap.index

# DROP NaN ...
df_cap = df_cap.dropna()

# start data analisis and signal extraction

X = df_cap.loc[:,:]
X_mean = X.mean(axis = 0).sum() 
X /= X_mean 

window_1 = 15  
window_2 = 30
         # Equal weighted
Z1 = X.rolling(window=window_1).mean()
Z2 = X.rolling(window=window_2).mean()

# remove NaN
X = X.iloc[window_2:,:]
Z1 = Z1.iloc[window_2:,:]
Z2 = Z2.iloc[window_2:,:]

date_from = '2011-01-01'
date_to = '2012-09-01' 

#-> Retrieve X, Z1 and Z2 only for the estimation window ...
X_window = X.loc[date_from:date_to, :]
Z1_window = Z1.loc[date_from:date_to, :]
Z2_window = Z2.loc[date_from:date_to, :]

params=[0.2,0.2,0.2,0.2,0.2]
lb = np.array([-0.01, 0.001, -5, 0.0001, 0.0001])
ub = np.array([5,   0.5,     5, 1,      1])
bounds = [(lb[0], ub[0]), (lb[1], ub[1]), (lb[2], ub[2]), (lb[3], ub[3]), (lb[4], ub[4])]
# PERFORM CALIBRATION
asset_id=1
res = minimize(mleMarket,params,args=(X_window,Z1_window,Z2_window,),bounds = bounds,
                       method='L-BFGS-B', options={'gtol': 1e-10, 'disp': True}, tol = 1e-9)
print(res)
# Creating an empty Dataframe with column names only
dfObj = pd.DataFrame(columns=['k', 'sigma', 'W0','W1','W2'],index=['S&P500'])
# Append rows in Empty Dataframe by adding dictionaries
dfObj.loc['S&P500']=[res.x[0],res.x[1],res.x[2],res.x[3],res.x[4]]
 
#print("Dataframe Contens ", dfObj, sep='\n')
dfObj<jupyter_output>fun: -532287.231423773
 hess_inv: <5x5 LbfgsInvHessProduct with dtype=float64>
      jac: array([-3.25962901, 35.01772881, -1.93249434, -1.68802217, -2.0256266 ])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 348
      nit: 40
   status: 0
  success: True
        x: array([0.02308758, 0.02946063, 0.02098286, 0.25364206, 0.25379852])
{'fraction_non_alphanumeric': 0.17851699094414059, 'fraction_numerical': 0.06330135389581279, 'mean_word_length': 3.5286236297198537, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '659117', 'license': 'MIT', 'n_tokens_mistral': 5010, 'n_tokens_neox': 4616, 'n_words': 905, 'path': 'IRL_market_model_DJI_calibration.ipynb', 'repo_name': 'Benn421/Coursera_Capstone', 'size': '11153'}