regression_3D.py

import numpy as np
import pandas as pd

from itertools import combinations
from stats import ols_tstats
from statsmodels.stats.stattools import durbin_watson
import matplotlib.pyplot as plt

def main(headlist):
  # Define experts names (Jackie Chan, Brad, Shahrukh)
  # expert_1 = "Shahrukh"
  # expert_2 = "Jackie Chan"

  # Load a specific column, e.g., column named "Sales"
  h1 = df[headlist[0]]  # Replace "Sales" with your column name
  h2 = df[headlist[1]]  # Replace "Sales" with your column name
  h3 = df[headlist[2]]  # Replace "Sales" with your column name
  h4 = df[headlist[3]]  # Replace "Sales" with your column name
  h5 = df[headlist[4]]  # Replace "Sales" with your column name
  h6 = df[headlist[5]]  # Replace "Sales" with your column name
  h7 = df[headlist[6]]  # Replace "Sales" with your column name
  h8 = df[headlist[7]]  # Replace "Sales" with your column name
  y = df[headlist[8]]

  n = len(h1)

  # ------------------------------------------------------------
  # 2. Design matrix with intercept
  # ------------------------------------------------------------
  H = np.column_stack([np.ones(n), h1, h2, h3, h4, h5, h6, h7, h8])  # shape (n,3)

  # Penalty matrix: intercept unpenalized
  L = np.diag([0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

  # ------------------------------------------------------------
  # 3. Ordinary Least Squares
  # ------------------------------------------------------------
  
  XtX_inv = np.linalg.inv(H.T @ H)
  theta_ols = XtX_inv @ (H.T @ y)
  yhat_ols = H @ theta_ols
  
  residuals = yhat_ols - y
  plt.figure(figsize=(6,4))
  plt.plot(range(n), residuals, color='purple', alpha=0.7, marker='o')
  plt.axhline(0, color='black', linestyle='--')
  plt.xlabel("Time")
  plt.ylabel("Residuals")
  plt.title("Residuals over Time")
  plt.grid(True)
  plt.tight_layout()
  plt.savefig("outputs/residuals_plot.png")
  plt.show()
  
  dw_stat = durbin_watson(residuals)
  print(f"Durbin-Watson statistic: {dw_stat:.4f}")
  
  # Plot residuals vs lagged residuals
  residuals_lagged = np.roll(residuals, 1)
  residuals_lagged[0] = 0  # Set first lagged value to 0 or np.nan

  plt.figure(figsize=(6,4))
  plt.scatter(residuals_lagged[1:], residuals[1:], color='teal', alpha=0.7)
  plt.xlabel("Residuals (lagged by 1 period)")
  plt.ylabel("Residuals")
  plt.title("Residuals vs Lagged Residuals")
  plt.grid(True)
  plt.tight_layout()
  plt.savefig("outputs/residuals_vs_lagged.png")
  plt.show()
  
  # Plot residuals vs h1 (PDI), h3 (PRICE), h5 (INVESTMENT), h6 (ADVERTISEMENT)
  features = [h1, h3, h5, h6]
  feature_names = [headlist[0], headlist[2], headlist[4], headlist[5]]

  plt.figure(figsize=(16, 4))
  for i, (feature, name) in enumerate(zip(features, feature_names)):
    plt.subplot(1, 4, i + 1)
    plt.scatter(feature, residuals, alpha=0.7, color='coral')
    plt.xlabel(name)
    plt.ylabel("Residuals")
    plt.title(f"Residuals vs {name}")
    plt.grid(True)
  plt.tight_layout()
  plt.savefig("outputs/residuals_vs_features.png")
  plt.show()
  
  plt.figure(figsize=(6,4))
  plt.scatter(yhat_ols, residuals, color='navy', alpha=0.7)
  plt.xlabel("Fitted Values (yhat_ols)")
  plt.ylabel("Residuals")
  plt.title("Residuals vs Fitted Values")
  plt.grid(True)
  plt.tight_layout()
  plt.savefig("outputs/residuals_vs_fitted.png")
  plt.show()
  
  plt.figure(figsize=(6,4))
  count, bins, ignored = plt.hist(residuals, bins=20, density=True, alpha=0.6, color='skyblue', edgecolor='black', label='Residuals Histogram')

  # Fit normal distribution
  mu, sigma = np.mean(residuals), np.std(residuals)
  x = np.linspace(bins[0], bins[-1], 100)
  plt.plot(x, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp(- (x - mu)**2 / (2 * sigma**2)), color='red', lw=2, label='Normal Distribution')

  plt.xlabel("Residuals")
  plt.ylabel("Density")
  plt.title("Histogram of Residuals with Normal Distribution")
  plt.legend()
  plt.grid(True)
  plt.tight_layout()
  plt.savefig("outputs/residuals_histogram.png")
  plt.show()
  
  sum_squared_error_ols = np.sum((yhat_ols - y)**2)

  print(f"OLS sum of square error: {sum_squared_error_ols}")

  # Display coefficients
  print([f"{header} - " for header in headlist])
  names = ["Intercept"] + [f"{header} - " for header in headlist]
  print("\nCoefficients:")
  for name, ols_c in zip(names, theta_ols):
      print(f"{name:10s}  OLS={ols_c:8.3f}")

  # ------------------------------------------------------------
  # 6. Plots
  # ------------------------------------------------------------

  # # (b) Predicted vs true scatter
  # plt.figure(figsize=(5,5))
  # plt.scatter(y, yhat_ols, c='blue', label=f'OLS, SE={np.round(sum_squared_error_ols, 5)}')
  # lims = [y.min()-1, y.max()+1]
  # plt.plot(lims, lims, 'k--')  # 45° reference
  # plt.xlabel("True y")
  # plt.ylabel("Predicted y")
  # plt.title(f"Predicted vs True, {expert_1} & {expert_2} & {expert_3}")
  # plt.legend()
  # plt.grid(True)
  # plt.savefig(f'pix_3D/predicted-vs-true_2D_{expert_1}-{expert_2}-{expert_3}.png')
  
  # plt.show()
  
  se_theta, t_stats, p_values = ols_tstats(y, H, theta_ols, XtX_inv)
  
  print("\nOLS computation statistics:")
  print("Coefficients:", theta_ols)
  print("Std Errors :", se_theta)
  print("t-stats    :", t_stats)
  print("p-values   :", p_values)
  
  # Build a DataFrame
  results_df = pd.DataFrame({
      "Header": ["intercept"] + headers[:-1],
      "theta_ols": theta_ols,
      "SE_Theta": se_theta,
      "T_Stats": t_stats,
      "P_Values": p_values
  })

  # Save to Excel
  results_df.to_excel("outputs/OLS_Results.xlsx", index=False)

  print("Results saved to OLS_Results.xlsx")

if __name__ == "__main__":
  
  # ------------------------------------------------------------
  # 1. Load data 
  # ------------------------------------------------------------

  # Load the Excel file
  file_path = "SALES.xlsx"
  df = pd.read_excel(file_path)  # Reads the first sheet by default

  # Display the first few rows to see the data
  print(df.head())

  headers = df.columns.tolist()
  # print(headers)
  
  main(headers)
  print('Execution done.')