API tools faq
paste
mirosh111000

Прикладна_економетрика_ЛР№2_Мірошниченко

mirosh111000
Oct 12th, 2025
146
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.49 KB | None | 0 0
raw download clone embed print report
  1. import math
  2. import numpy as np
  3. import pandas as pd
  4. import matplotlib.pyplot as plt
  5.  
  6. pd.set_option('display.width', 120)
  7. pd.set_option('display.max_columns', 20)
  8. pd.set_option('display.precision', 4)
  9. from scipy.stats import t as t_dist, chi2 as chi2_dist
  10.  
  11. x = np.array([0.3, 0.6, 1.1, 0.2, 0.8, 1.6, 1.5, 0.9, 0.9, 0.8], dtype=float)
  12. y = np.array([2.8, 4.3, 6.3, 2.9, 5.3, 8.3, 8.0, 5.7, 5.4, 5.2], dtype=float)
  13. x_new = 2.0
  14.  
  15. i = np.arange(1, len(y) + 1, dtype=int)
  16. T1 = pd.DataFrame({
  17.     'i': i,
  18.     'x_i': x,
  19.     'y_i': y,
  20.     'x_i^2': x ** 2,
  21.     'x_i*y_i': x * y,
  22.     'y_i^2': y ** 2,
  23. })
  24. T1_sum = pd.DataFrame({
  25.     'i': ['Σ'],
  26.     'x_i': [x.sum()],
  27.     'y_i': [y.sum()],
  28.     'x_i^2': [(x ** 2).sum()],
  29.     'x_i*y_i': [(x * y).sum()],
  30.     'y_i^2': [(y ** 2).sum()],
  31. })
  32. T1_print = pd.concat([T1.round(6), T1_sum], ignore_index=True)
  33. print('\n=== Таблиця 1: ===')
  34. print(T1_print.to_string(index=False))
  35.  
  36. x_bar = x.mean()
  37. y_bar = y.mean()
  38. Sxx = float(np.sum((x - x_bar) ** 2))
  39. Sxy = float(np.sum((x - x_bar) * (y - y_bar)))
  40.  
  41. a1 = Sxy / Sxx
  42. a0 = y_bar - a1 * x_bar
  43.  
  44. y_hat = a0 + a1 * x
  45. u = y - y_hat
  46. u2 = u ** 2
  47.  
  48. T2 = pd.DataFrame({
  49.     'i': i,
  50.     'y_i': y,
  51.     'ŷ_i': y_hat,
  52.     'u_i = y_i - ŷ_i': u,
  53.     'u_i^2': u2,
  54. })
  55. T2_sum = pd.DataFrame({
  56.     'i': ['Σ'],
  57.     'y_i': [''],
  58.     'ŷ_i': [''],
  59.     'u_i = y_i - ŷ_i': [u.sum()],
  60.     'u_i^2': [u2.sum()],
  61. })
  62. T2_print = pd.concat([T2.round(6), T2_sum], ignore_index=True)
  63. print('\n=== Таблиця 2: ===')
  64. print(T2_print.to_string(index=False))
  65.  
  66. RSS = float(u2.sum())
  67. n = len(x)
  68. k = n - 2
  69. s2 = RSS / k
  70. s = math.sqrt(s2)
  71.  
  72. alpha_90 = 0.10
  73. chi2_low = float(chi2_dist.ppf(1 - alpha_90/2, k))
  74. chi2_high = float(chi2_dist.ppf(alpha_90/2, k))
  75.  
  76. sigma2_CI_low = (k * s2) / chi2_low
  77. sigma2_CI_high = (k * s2) / chi2_high
  78.  
  79. sigma2_df = pd.DataFrame({
  80.     's2 (оцінка σ²)': [s2],
  81.     's = sqrt(s2)': [s],
  82.     'df (k)': [k],
  83.     'CI90_low': [sigma2_CI_low],
  84.     'CI90_high': [sigma2_CI_high],
  85. }).round(6)
  86. print('\n=== Дисперсія збурень σ²: точкова та 90% ДІ ===')
  87. print(sigma2_df.to_string(index=False))
  88.  
  89. SE_a1 = s / math.sqrt(Sxx)
  90. SE_a0 = s * math.sqrt(1/n + (x_bar**2)/Sxx)
  91. alpha_95 = 0.05
  92. t_crit = float(t_dist.ppf(1 - alpha_95/2, k))
  93.  
  94. t_a0 = a0 / SE_a0
  95. t_a1 = a1 / SE_a1
  96. p_a0 = 2 * (1 - t_dist.cdf(abs(t_a0), df=k))
  97. p_a1 = 2 * (1 - t_dist.cdf(abs(t_a1), df=k))
  98.  
  99. a0_CI_low = a0 - t_crit * SE_a0
  100. a0_CI_high = a0 + t_crit * SE_a0
  101. a1_CI_low = a1 - t_crit * SE_a1
  102. a1_CI_high = a1 + t_crit * SE_a1
  103.  
  104. coef_df = pd.DataFrame({
  105.     'coef': ['a0', 'a1'],
  106.     'estimate': [a0, a1],
  107.     'SE': [SE_a0, SE_a1],
  108.     't': [t_a0, t_a1],
  109.     'p_value(two‑tailed)': [p_a0, p_a1],
  110.     'CI95_low': [a0_CI_low, a1_CI_low],
  111.     'CI95_high': [a0_CI_high, a1_CI_high],
  112. }).round(4)
  113. print('\n=== Коефіцієнти моделі: оцінки, SE, t, p, 95% ДІ ===')
  114. print(coef_df.to_string(index=False))
  115.  
  116. TSS = float(np.sum((y - y_bar) ** 2))
  117. R2 = 1 - RSS / TSS
  118. r = math.copysign(math.sqrt(R2), a1)
  119. MAPE = float(np.mean(np.abs((y_hat - y) / y) * 100))
  120. MPE = float(np.mean(((y_hat - y) / y) * 100))
  121. metrics_df = pd.DataFrame({
  122.     'R2': [R2],
  123.     'r = sign(a1)*sqrt(R2)': [r],
  124.     'MAPE (%)': [MAPE],
  125.     'MPE (%)': [MPE],
  126. }).round(4)
  127. print('\n=== Метрики якості моделі (R², r, MAPE, MPE) ===')
  128. print(metrics_df.to_string(index=False))
  129.  
  130. S_yhat_i = s * np.sqrt(1/n + (x - x_bar) ** 2 / (Sxx))
  131. upper_i = y_hat + t_crit * S_yhat_i
  132. lower_i = y_hat - t_crit * S_yhat_i
  133. order = np.argsort(x)
  134.  
  135. plt.figure(figsize=(7, 5))
  136. plt.scatter(x, y, label='Дані', zorder=3)
  137. plt.plot(x[order], y_hat[order], label='OLS лінія', linewidth=2, zorder=4)
  138. plt.plot(x[order], upper_i[order], linestyle='--', label='Верхня межа 95%', zorder=2)
  139. plt.plot(x[order], lower_i[order], linestyle='--', label='Нижня межа 95%', zorder=2)
  140. plt.xlabel('x')
  141. plt.ylabel('y')
  142. plt.title('Дані, OLS та 95% інтервал довіри')
  143. plt.legend()
  144. plt.grid(alpha=0.25)
  145. plt.tight_layout()
  146. plt.show()
  147.  
  148. SE_pred = s * math.sqrt(1 + 1/n + (x_new - x_bar) ** 2 / Sxx)
  149. y_new_hat = a0 + a1 * x_new
  150. PI_low = y_new_hat - t_crit * SE_pred
  151. PI_high = y_new_hat + t_crit * SE_pred
  152.  
  153. pred_df = pd.DataFrame({
  154.     'x_new': [x_new],
  155.     'ŷ(x_new)': [y_new_hat],
  156.     'SE_pred': [SE_pred],
  157.     'PI95_low': [PI_low],
  158.     'PI95_high': [PI_high],
  159. }).round(4)
  160. print('\n=== Прогноз у x_new та 95% ДІ прогнозу (PI) ===')
  161. print(pred_df.to_string(index=False))
  162.  
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!