<jupyter_start><jupyter_text>Create weighted average of RMSE:<jupyter_code>def weighted_avg (RMSEs):
#     print(RMSEs)
    sums = 0
    divide = 1
    #  Sum of all RMSE
    for i in range(1,len(RMSEs)+1):
        name = 'M{}'.format(str(i))
        sums  = sums + RMSEs[name] #sum1
        print(sums)
    sums2 = 0
    
    for i in range(1, len(RMSEs)+1):
        #  Devide total sum by individual RMSE Value
        name = 'M{}'.format(str(i))
        RMSEs[name] = sums/RMSEs[name] #n1
        #  Add all processed RMSE 
        sums2 = sums2 + RMSEs[name]
        
    #  Again We devide RMSE by total sum to get the ratio 
    for i in range(1, len(RMSEs)+1):
        name = 'M{}'.format(str(i))
        RMSEs[name] = RMSEs[name]/sums2
    
    print(RMSEs)
# {'df_ar':1.29, 'df_ema':4.24, 'df_macd':4.31, 'df_poly':19.6, 'df_LR':20.8}
dict_algo = {'M1':1.29, 'M2':4.24, 'M3':4.31, 'M4':19.6, 'M5':20.8}
weighted_avg(dict_algo)
df_ar['Predict_AutoReg'] = 0.5775723622356099 * df_ar['Predict_AutoReg']
df_ema['Predict_EMA'] = 0.17572366681224924 * df_ema['Predict_EMA']
df_macd['Predict_MACD'] = 0.17286968614476494 * df_macd['Predict_MACD']
df_poly['Predictions'] = 0.03801369118795596 * df_poly['Predictions']
df_LR['Preds'] = 0.03582059361942004 * df_LR['Preds']
y_ped = pd.DataFrame()
y_ped['close'] = df_ar['Predict_AutoReg'] + df_ema['Predict_EMA'] + df_macd['Predict_MACD'] + df_poly['Predictions'] + df_LR['Preds']
y_ped['close']
data = y_ped['close'].dropna()
data_to_use = data.values
data_to_use<jupyter_output><empty_output><jupyter_text>The idea behind StandardScaler is that it will transform your data such that its distribution will have a mean value 0 and standard deviation of 1.<jupyter_code>scaler = StandardScaler()
scaled_data = scaler.fit_transform(data_to_use.reshape(-1, 1))
scaled_data
# z = (x - u) / s
# scaled_data = (((data_to_use)-data_to_use.mean()) / data_to_use.std().reshape(-1,1)).reshape(-1, 1)
# scaled_data
plt.figure(figsize=(16,7), frameon=False, facecolor='brown', edgecolor='blue')
plt.title('Weighted Apple stocks from September 2018 to June 2019')
plt.xlabel('Days')
plt.ylabel('Scaled value of stocks')
plt.plot(scaled_data, label='Stocks data')
plt.legend()
plt.show()
def lstm_tf(scaledData):
    
#     Taking 5 X 5 size window each time
    def window_data(data, window_size):
        
        X = []
        y = []
        
        i = 0
        while (i + window_size) <= len(data) - 1:
            X.append(data[i:i+window_size])
            y.append(data[i+window_size])
            
            i += 1
        assert len(X) ==  len(y)
        return X, y
    
    X, y = window_data(scaled_data, 5)
#     print("X, y :", X, y)
    X_train  = np.array(X[:150])
    y_train = np.array(y[:150])
    
    X_test = np.array(X[150:])
    y_test = np.array(y[150:])
    print(" ----------------------- Splitting Data Into Train Test ----------------------- ")
    print("LSTM Model Information:")
    print("X_train size: {}".format(X_train.shape))
    print("y_train size: {}".format(y_train.shape))
    print("X_test size: {}".format(X_test.shape))
    print("y_test size: {}".format(y_test.shape))
    print("-------------------------------------------------------------------------------- ")

    #Hyperparameters used in the LSTM network
    #how many windows of data we are passing at once
    batch_size = 5
    #how big window_size is (Or How many days do we consider to predict next point in the sequence)
    window_size = 5     
    #How many units do we use in LSTM cell
    hidden_layer = 256 
    #To prevent exploding gradient, we use clipper to clip gradients below -margin or above this margin
    clip_margin = 4 
    learning_rate = 0.001 
    epochs = 150
    
    
    inputs = tf.placeholder(tf.float32, [batch_size, window_size, 1])
    targets = tf.placeholder(tf.float32, [batch_size, 1])
    
    # LSTM weights
    # Weights for the input gate
        # Outputs random values from a truncated normal distribution    
    weights_input_gate = tf.Variable(tf.truncated_normal([1, hidden_layer], stddev=0.05))
    weights_input_hidden = tf.Variable(tf.truncated_normal([hidden_layer, hidden_layer], stddev=0.05))
    bias_input = tf.Variable(tf.zeros([hidden_layer]))
    
    # weights for the forgot gate
    weights_forget_gate = tf.Variable(tf.truncated_normal([1, hidden_layer], stddev=0.05))
    weights_forget_hidden = tf.Variable(tf.truncated_normal([hidden_layer, hidden_layer], stddev=0.05))
    bias_forget = tf.Variable(tf.zeros([hidden_layer]))
    
    # weights for the output gate
    weights_output_gate = tf.Variable(tf.truncated_normal([1, hidden_layer], stddev=0.05))
    weights_output_hidden = tf.Variable(tf.truncated_normal([hidden_layer, hidden_layer], stddev=0.05))
    bias_output = tf.Variable(tf.zeros([hidden_layer]))
    
    # weights for the memory cell
    weights_memory_cell = tf.Variable(tf.truncated_normal([1, hidden_layer], stddev=0.05))
    weights_memory_cell_hidden = tf.Variable(tf.truncated_normal([hidden_layer, hidden_layer], stddev=0.05))
    bias_memory_cell = tf.Variable(tf.zeros([hidden_layer]))
    
    ## Output layer weigts
    weights_output = tf.Variable(tf.truncated_normal([hidden_layer, 1], stddev=0.05))
    bias_output_layer = tf.Variable(tf.zeros([1]))
    
    
    def LSTM_cell(input, output, state):
        
        input_gate = tf.sigmoid(tf.matmul(input, weights_input_gate) + tf.matmul(output, weights_input_hidden) + bias_input)
        
        forget_gate = tf.sigmoid(tf.matmul(input, weights_forget_gate) + tf.matmul(output, weights_forget_hidden) + bias_forget)
        
        output_gate = tf.sigmoid(tf.matmul(input, weights_output_gate) + tf.matmul(output, weights_output_hidden) + bias_output)
        
        memory_cell = tf.tanh(tf.matmul(input, weights_memory_cell) + tf.matmul(output, weights_memory_cell_hidden) + bias_memory_cell)
        
        state = state * forget_gate + input_gate * memory_cell
        
        output = output_gate * tf.tanh(state)
        return state, output
    
    
    outputs = []
    for i in range(batch_size): #Iterates through every window in the batch
        #for each batch I am creating batch_state as all zeros and output for that window which is all zeros at the beginning as well.
        batch_state = np.zeros([1, hidden_layer], dtype=np.float32) 
        batch_output = np.zeros([1, hidden_layer], dtype=np.float32)
        #for each point in the window we are feeding that into LSTM to get next output
        for ii in range(window_size):
            batch_state, batch_output = LSTM_cell(tf.reshape(inputs[i][ii], (-1, 1)), batch_state, batch_output)
        #last output is conisdered and used to get a prediction
        outputs.append(tf.matmul(batch_output, weights_output) + bias_output_layer)
        
        
    outputs
    
    losses = []
    
    for i in range(len(outputs)):
        losses.append(tf.losses.mean_squared_error(tf.reshape(targets[i], (-1, 1)), outputs[i]))
        
    loss = tf.reduce_mean(losses)
    
    gradients = tf.gradients(loss, tf.trainable_variables())
    clipped, _ = tf.clip_by_global_norm(gradients, clip_margin)
    optimizer = tf.train.AdamOptimizer(learning_rate)
    trained_optimizer = optimizer.apply_gradients(zip(gradients, tf.trainable_variables()))
    
    session = tf.Session()
    session.run(tf.global_variables_initializer())
    
    for i in range(epochs):
        traind_scores = []
        ii = 0
        epoch_loss = []
        while(ii + batch_size) <= len(X_train):
            X_batch = X_train[ii:ii+batch_size]
            y_batch = y_train[ii:ii+batch_size]
            
            o, c, _ = session.run([outputs, loss, trained_optimizer], feed_dict={inputs:X_batch, targets:y_batch})
            
            epoch_loss.append(c)
            traind_scores.append(o)
            ii += batch_size
        if (i % 10) == 0:
            print('Epoch {}/{}'.format(i, epochs), ' Current loss: {}'.format(np.mean(epoch_loss)))
    
    
    
    sup =[]
    for i in range(len(traind_scores)):
        for j in range(len(traind_scores[i])):
            sup.append(traind_scores[i][j][0])
            
    print("len(X_test): ",len(X_test))      
    tests = []
    i = 0
    while i+batch_size <= len(X_test):
        
        o = session.run([outputs], feed_dict={inputs:X_test[i:i+batch_size]})
        i += batch_size
        tests.append(o)
        
    print('length of text: ',len(tests))
    tests_new = []
    for i in range(len(tests)):
        for j in range(len(tests[i][0])):
            tests_new.append(tests[i][0][j])
            
    test_results = []
    for i in range(77):
#         print("i: ",i)
        if i >= 40:
            test_results.append(tests_new[i-40])
        else:
            test_results.append(None)
    
    plt.figure(figsize=(16, 7))
    plt.plot(scaled_data, label='Original data')
    plt.plot(sup, label='Training data')
    plt.plot(range(113,190),test_results, label='Testing data')
    plt.legend()
    plt.show()
    
    inv_ytrain = scaler.inverse_transform(np.array(y_test).reshape(-1, 1))
    inv_test = scaler.inverse_transform(np.array(tests_new).reshape(-1, 1))
#     print("size inv_test: ", inv_test.shape)
#     print("size inv_ytrain: ", inv_ytrain[0:40].shape)
    
    residuals = inv_test - inv_ytrain[0:40]
    rmse = np.sqrt(np.mean(residuals ** 2))
    R2 =1-((residuals.transpose()@residuals)/(np.var(inv_ytrain)*len(sup)))
    print('RMSE is: ', rmse)   
    session.close()
    
    #indata = scaler.inverse_transform(test_results)
    
    return R2, inv_test, inv_ytrain
Rsq, pred_values, test_values = lstm_tf(scaled_data)
# Zooming in on the prediction and actual
fig , ax = plt.subplots(figsize =(15,6))
ax.plot(test_values, color='red', label='True Closing Price')
plt.plot(pred_values, color='green', label='Predicted Closing Price' )<jupyter_output><empty_output><jupyter_text>Calculating Residuals<jupyter_code>residual_errors = (test_values[:40] - pred_values)<jupyter_output><empty_output><jupyter_text>Plotting error bar plots<jupyter_code>from matplotlib.pyplot import figure
figure(num=None, figsize=(14, 8), dpi=80, facecolor='w', edgecolor='red')
plt.errorbar(np.arange(len(test_values[:40])),test_values[:40],yerr=residual_errors,uplims=True, lolims=True)<jupyter_output>C:\Users\kaila\Anaconda3\envs\tensorflow_gpuenv\lib\site-packages\ipykernel_launcher.py:3: MatplotlibDeprecationWarning: Support for passing a (n, 1)-shaped error array to errorbar() is deprecated since Matplotlib 3.1 and will be removed in 3.3; pass a 1D array instead.
  This is separate from the ipykernel package so we can avoid doing imports until<jupyter_text>Plotting QQ plot to see if the errors are normally distributed.<jupyter_code>from scipy import stats
plt.figure(figsize=(14, 8))
res = stats.probplot(residual_errors.flatten(), plot=plt)
plt.show()<jupyter_output><empty_output><jupyter_text>- From the QQ plot we can see that the errors have a curve and does nothave a linear fit.<jupyter_code>from sklearn.ensemble.forest import RandomForestRegressor
def modelfit(alg, dtrain, predictors, useTrainCV=True, cv_folds=5, early_stopping_rounds=50):    
    if useTrainCV:
        xgb_param = alg.get_xgb_params()
        xgtrain = xgb.DMatrix(dtrain[predictors].values, label=dtrain[target].values)
        cvresult = xgb.cv(xgb_param, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds,
            metrics='rmse', early_stopping_rounds=early_stopping_rounds, show_stdv=False, verbose_eval=False)
        alg.set_params(n_estimators=cvresult.shape[0])
            
    #Fit the algorithm on the data
    alg.fit(dtrain[predictors], dtrain['median_apt_value'],eval_metric='rmse')
    #Predict training set:
    dtrain_predictions = alg.predict(dtrain[predictors])
    print("Cross validation results: ")
    print("Train RMSE Min: ",round(cvresult['train-rmse-mean'].min(),4), "      Train RMSE Max: ",round(cvresult['train-rmse-mean'].max(),4))
    print("Test RMSE Min: ",round(cvresult['test-rmse-mean'].min(),4), "      Test RMSE Max: ",round(cvresult['test-rmse-mean'].max(),4))
    
    print ("R2 Score of model : %.4g" % metrics.r2_score(df_pricing['median_apt_value'].values, dtrain_predictions))
    
target="median_apt_value"
predictors = [x for x in df_pricing.columns if x not in [target]]

for mx_dpt in range(2,12,2):
    print('\n------------------- Prediction at max_depth= ',mx_dpt,'-------------------')
    xgb2 = xgboost.XGBRegressor(
        max_depth=mx_dpt,
        learning_rate =0.07,
        n_estimators=500,
        objective= 'reg:squarederror',
        min_child_weight=1.2,
        gamma=0.05,
        subsample=0.7,
        colsample_bytree=0.5,
        # nthread=11,
        n_jobs=11,
        # scale_pos_weight=1,
        # reg_alpha=0.1,
         seed=45)
    modelfit(xgb2, df_pricing, predictors)
    
# build our RF model
RF_Model = RandomForestRegressor(n_estimators=100, max_features=1, oob_score=True)
 
# let's get the labels and features in order to run our 
# model fitting
labels = y_train#[:, None]
features = X_train[:, None]
 
# Fit the RF model with features and labels.
rgr=RF_Model.fit(features, labels)
 
# Now that we've run our models and fit it, let's create
# dataframes to look at the results
X_test_predict=pd.DataFrame(
    rgr.predict(X_test[:, None])).rename(
    columns={0:'predicted_price'}).set_index('predicted_price')
X_train_predict=pd.DataFrame(
    rgr.predict(X_train[:, None])).rename(
    columns={0:'predicted_price'}).set_index('predicted_price')
 
# combine the training and testing dataframes to visualize
# and compare.
RF_predict = X_train_predict.append(X_test_predict)<jupyter_output><empty_output><jupyter_text>Generating Additional data for Multivariat<jupyter_code>from pandas import DataFrame
from pandas import concat
 
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
    n_vars = 1 if type(data) is list else data.shape[1]
    df = DataFrame(data)
    cols, names = list(), list()
    # input sequence (t-n, ... t-1)
    for i in range(n_in, 0, -1):
        cols.append(df.shift(i))
        names += [('close%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
    # forecast sequence (t, t+1, ... t+n)
    for i in range(0, n_out):
        cols.append(df.shift(-i))
        if i == 0:
            names += [('close%d(t)' % (j+1)) for j in range(n_vars)]
        else:
            names += [('close%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
    # put it all together
    agg = concat(cols, axis=1)
    agg.columns = names
    # drop rows with NaN values
    if dropnan:
        agg.dropna(inplace=True)
    return agg

values = scaled_data
data = series_to_supervised(values, 3)
print(data.tail(20))<jupyter_output>close1(t-3)  close1(t-2)  close1(t-1)  close1(t)
858     0.653056     0.677579     0.591482   0.544302
859     0.677579     0.591482     0.544302   0.713297
860     0.591482     0.544302     0.713297   0.790597
861     0.544302     0.713297     0.790597   0.862034
862     0.713297     0.790597     0.862034   0.993445
863     0.790597     0.862034     0.993445   1.058218
864     0.862034     0.993445     1.058218   1.117659
865     0.993445     1.058218     1.117659   1.101133
866     1.058218     1.117659     1.101133   1.100066
867     1.117659     1.101133     1.100066   1.062483
868     1.101133     1.100066     1.062483   1.093136
869     1.100066     1.062483     1.093136   1.214684
870     1.062483     1.093136     1.214684   1.199224
871     1.093136     1.214684     1.199224   1.241607
872     1.214684     1.199224     1.241607   1.223481
873     1.199224     1.241607     1.223481   1.218150
874     1.241607     1.223481     1.218150   1.137917
875     1.223481     1.21815[...]<jupyter_text>Trading Strategy  Bollinger_band- A Bollinger Band is a technical analysis tool defined by a set of lines plotted two standard deviations (positively and negatively) away from a simple moving average (SMA)- Bollinger Bands are a technical analysis tool developed by <NAME>.- There are three lines that compose Bollinger Bands: A simple moving average (middle band) and an upper and lower band.- The upper and lower bands are typically 2 standard deviations +/- from a 20-day simple moving average, but can be modified.<jupyter_code>#read csv file, use date as index and read close as a column
df = pd.read_csv('AAPL.csv'.format('AAPL'), index_col='Date', parse_dates=True, usecols=['Date','Close'])
print(df)
df = df.rename(columns={'Close': 'AAPL STOCK'})

# calculate Simple Moving Average with 20 days window
sma = df.rolling(20, win_type='triang').mean()

# calculate the standard deviation
rstd = df.rolling(20).std()

upper_band = sma + 2 * rstd
upper_band = upper_band.rename(columns={'AAPL STOCK': 'upper'})
df['Upper Band'] = upper_band.upper
# print(upper_band)

lower_band = sma - 2 * rstd
lower_band = lower_band.rename(columns={'AAPL STOCK': 'lower'})
df['Lower Band'] = lower_band.lower
# print(lower_band)


ax = df.plot(title='{} Price and Bollinger Band'.format('AAPL STOCK'), figsize=(17,8))
ax.fill_between(df.index, df['Lower Band'], df['Upper Band'], color='#ADCCFF', alpha='0.4')
ax.set_xlabel('Date')
ax.set_ylabel('SMA and BB')
ax.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>- In the chart depicted above, Bollinger Bands bracket the 20-day SMA of the stock with an upper and lower band along with the daily movements of the stock's price. Because standard deviation is a measure of volatility, when the markets become more volatile the bands widen; during less volatile periods, the bands contract.  Take weighted Average of RMSE Values: In statistics, inverse-variance weighting is a method of aggregating two or more random variables to minimize the variance of the weighted average. Each random variable is weighted in inverse proportion to its variance.Given a sequence of independent observations yi with variances σi2, the inverse-variance weighted average is given by[1]![image.png](attachment:image.png) The inverse-variance weighted average has the least variance among all weighted averages, which can be calculated as![image.png](attachment:image.png)<jupyter_code>rmse_vals = np.array([1.1,6.7,6.8,7.2,11.1,3.4])
mu = rmse_vals.mean()
print('mu: ',mu)
sigma2 = rmse_vals.var()
print('sigma square: ',sigma2)
normalised_value = (rmse_vals-mu)/sigma2
normalised_value
print(sum(rmse_vals))
wt_vals2 = []
new_vals = 0.0
for i,val in enumerate(wt_vals):
#     print(val)
    new_vals = 0.0
    new_vals = 1 / (1/sigma2)
    wt_vals2.append(new_vals)
    
print("wsvsv",wt_vals2)
wt_vals = (sum(rmse_vals)/sigma2) / (1/sigma2)
wt_vals
rmse_vals = [1.1,6.7,6.8,7.2,11.1,3.4]
wt_vals = []
for i,val in enumerate(rmse_vals):
#     print(val)
    new_vals = 0.0
    new_vals = (val/sigma2) / (1/sigma2)
    wt_vals.append(new_vals)
    
print("wsvsv",wt_vals)<jupyter_output>wsvsv [1.1, 6.699999999999999, 6.8, 7.2, 11.1, 3.4]<jupyter_text>Calculate percentage of RMSE of individual model:<jupyter_code>rmse_vals = np.array([1.1,6.7,6.8,7.2,11.1,3.4])
wt_perc = []
total_rmse = sum(rmse_vals)

for i,val in enumerate(sorted(rmse_vals)):
    per_val = 0.0
    per_val = total_rmse / 
    wt_perc.append(round(per_val,4))
    
print("Percentage of RMSE: ",wt_perc)<jupyter_output>Percentage of RMSE:  [0.0303, 0.0937, 0.1846, 0.1873, 0.1983, 0.3058]<jupyter_text>Checking if total value is adding upto 1<jupyter_code>sum(wt_perc)<jupyter_output><empty_output><jupyter_text>Assigning highest value to model with low RMSE:<jupyter_code># revesing the array will simply give the highes value to 1st
wt_inv_perc = wt_perc[::-1]
print("Reversed weighted average: ",wt_inv_perc)
sum(wt_inv_perc)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1332119556971628, 'fraction_numerical': 0.05345673393010671, 'mean_word_length': 3.399110122358176, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '267081', 'license': 'MIT', 'n_tokens_mistral': 7642, 'n_tokens_neox': 6956, 'n_words': 1590, 'path': 'APPL_LSTM.ipynb', 'repo_name': 'KailashDN/Statistical_modeling_AAPL_stock_prediction', 'size': '19781'}