# Energy price non-linear regression # solve for oil sales price (outcome) # using 3 predictors of WTI Oil Price, # Henry Hub Price and MB Propane Spot Price import numpy as np from scipy.optimize import minimize import pandas as pd import numpy as np import matplotlib.pyplot as plt # data file from URL address data = 'https://apmonitor.com/me575/uploads/Main/oil_data.txt' df = pd.read_csv(data) xm1 = np.array(df["WTI_PRICE"]) # WTI Oil Price xm2 = np.array(df["HH_PRICE"]) # Henry Hub Gas Price xm3 = np.array(df["NGL_PRICE"]) # MB Propane Spot Price ym = np.array(df["BEST_PRICE"]) # oil sales price received (outcome) # calculate y def calc_y(x): a = x[0] b = x[1] c = x[2] d = x[3] #y = a * xm1 + b # linear regression y = a * ( xm1 ** b ) * ( xm2 ** c ) * ( xm3 ** d ) return y # define objective def objective(x): # calculate y y = calc_y(x) # calculate objective obj = 0.0 for i in range(len(ym)): obj = obj + ((y[i]-ym[i])/ym[i])**2 # return result return obj # initial guesses x0 = np.zeros(4) x0[0] = 0.0 # a x0[1] = 0.0 # b x0[2] = 0.0 # c x0[3] = 0.0 # d # show initial objective print('Initial Objective: ' + str(objective(x0))) # optimize # bounds on variables my_bnds = (-100.0, 100.0) bnds = (my_bnds, my_bnds, my_bnds, my_bnds) solution = minimize(objective, x0, method='SLSQP', bounds=bnds) x = solution.x y = calc_y(x) # show final objective cObjective = 'Final Objective: ' + str(objective(x)) print(cObjective) # print solution print('Solution') cA = 'A = ' + str(x[0]) print(cA) cB = 'B = ' + str(x[1]) print(cB) cC = 'C = ' + str(x[2]) print(cC) cD = 'D = ' + str(x[3]) print(cD) cFormula = "Formula is : " + "\n" \ + "A * WTI^B * HH^C * PROPANE^D" cLegend = cFormula + "\n" + cA + "\n" + cB + "\n" \ + cC + "\n" + cD + "\n" + cObjective #ym measured outcome #y predicted outcome from scipy import stats slope, intercept, r_value, p_value, std_err = stats.linregress(ym,y) r2 = r_value**2 cR2 = "R^2 correlation = " + str(r_value**2) print(cR2) # plot solution plt.figure(1) plt.title('Actual (YM) versus Predicted (Y) Outcomes For Non-Linear Regression') plt.plot(ym,y,'o') plt.xlabel('Measured Outcome (YM)') plt.ylabel('Predicted Outcome (Y)') plt.legend([cLegend]) plt.grid(True) plt.show()