evgeniya_polyntseva

smad 4

Dec 1st, 2020
501
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. from statistics import mean
  3. from scipy.stats import chi2
  4. import pandas as pd
  5. from scipy.stats import f
  6. import matplotlib.pyplot as plt
  7.  
  8.  
  9. def func(x1, x2):
  10.     return 1, 1 / x1, x1, x2 ** 2, x2
  11.  
  12.  
  13. def func_u(x1, x2, theta):
  14.     return theta[0] * func(x1, x2)[0] + theta[1] * func(x1, x2)[1] + theta[2] * func(x1, x2)[2] + theta[3] * \
  15.            func(x1, x2)[3] + \
  16.            theta[4] * func(x1, x2)[4]
  17.  
  18.  
  19. def point_generation(a, b, n):
  20.     rng = np.random.default_rng()
  21.     p1 = (b - a) * rng.random(n) + a
  22.     p2 = (b - a) * rng.random(n) + a
  23.     return p1, p2
  24.  
  25.  
  26. def simulation(u, n, m, p, x1, x2):
  27.     omega = np.dot((u - mean(u)), (u - mean(u)).T) / (n - 1)
  28.     disp = np.math.sqrt(p * omega)
  29.     y = u + np.random.normal(0, disp, n)
  30.     # theta estimate
  31.     x = np.array([np.ones(n), func(x1, x2)[1], func(x1, x2)[2], func(x1, x2)[3], func(x1, x2)[4]])
  32.     theta_est = np.linalg.inv(x.dot(x.T)).dot(x).dot(y)
  33.     # sigma estimate
  34.     u_est = func_u(x1, x2, theta_est)
  35.     e = y - u_est
  36.     sigma_est = np.dot(e, e.T) / (n - m)
  37.     return y, x, theta_est, sigma_est
  38.  
  39.  
  40. # Генерация гетероскедатичности
  41. def heteroscedaticity(u, n):
  42.     y_gen = []
  43.     for i in range(n):
  44.         y_gen.append(u[i] + np.random.normal(0, np.exp(abs(u[i]))))
  45.     return y_gen
  46.  
  47.  
  48. # Функция для генерации вектора известных переменных
  49. def known_var(var, n):
  50.     v = []
  51.     for i in range(n):
  52.         v.append(float(np.exp(abs(var[i]))))
  53.     v = np.array(v)
  54.     return v
  55.  
  56.  
  57. # 2. Полученные данные проверить по тестам на наличие гетероскедастичности.
  58. # 1) Тест Бреуша-Пагана (не работает, вид заебись)
  59. def Breusha_Pagana(u, n):
  60.     y_t = np.array(heteroscedaticity(u, n))
  61.     # МНК и ОМП
  62.     e_t = y_t - u
  63.     sigma_tilda = sum((e_t ** 2) / n)
  64.     c_t = e_t ** 2 / sigma_tilda
  65.     # Построение регрессии и вычисление ESS
  66.     z = np.array([np.ones(n), known_var(y_t, n)])
  67.     alpha = (np.linalg.inv(z.dot(z.T)).dot(z)).dot(c_t.T)
  68.     c_i = z.T * alpha
  69.     c_mean = c_i.mean()
  70.     ESS = 0
  71.     for i in range(n):
  72.         ESS += (c_i[i] - c_mean) ** 2
  73.     a = ESS / 2
  74.     if a < chi2.ppf(0.95, 2):
  75.         print('Гипотеза о гомоскедастичности принимается.')
  76.     else:
  77.         print('Гипотеза о гомоскедастичности отвергается.')
  78.  
  79.  
  80. # 2) Тест Голдфельда-Квандтона (работает, нужно привести код в человеческий вид)
  81. def Goldfeld_Quandton(u, n, m, p, x1, x2):
  82.     k = m
  83.     y, x, theta_est, sigma_est = simulation(u, n, m, p, x1, x2)
  84.     s = known_var(y, n)
  85.     values = pd.DataFrame({'x1': x1, 'x2': x2, 'y': y, 'u': u, 'exp(|y|)': s})
  86.     # Сортировка по возрастанию
  87.     values = values.sort_values('exp(|y|)')
  88.     values = values.reset_index(drop=True)
  89.     pd.set_option("display.max_rows", None, "display.max_columns", None)
  90.     nc = int(n / 3)
  91.     # Оценка по первым (n-nc)/2 наблюдений
  92.     newTable = values
  93.     X = []
  94.     for i in range(len(theta_est)):
  95.         X.append([])
  96.         for j in range(nc):
  97.             X[i].append(func(newTable['x1'][j], newTable['x2'][j])[i])
  98.     X = np.matrix(X)
  99.     Q_T = np.linalg.inv(X.dot(X.T))
  100.     Q_T = Q_T.dot(X)
  101.     Q_T = Q_T.dot(np.matrix(newTable['y'][0:nc]).T)
  102.     RSS1 = (np.array(newTable['y'][0:nc]) - X.T.dot(Q_T).T)
  103.     RSS1 = float(RSS1.dot(RSS1.T))
  104.     # Оценка по последним (n-nc)/2 наблюдений
  105.     X = []
  106.     for i in range(len(theta_est)):
  107.         X.append([])
  108.         for j in range(nc):
  109.             X[i].append(func(newTable['x1'][j + 2 * nc], newTable['x2'][j + 2 * nc])[i])
  110.     X = np.matrix(X)
  111.     Q_T = np.linalg.inv(X.dot(X.T))
  112.     Q_T = Q_T.dot(X)
  113.     newTable = newTable.reset_index(drop=True)
  114.     Q_T = Q_T.dot(np.matrix(newTable.y[2 * nc:]).T)
  115.     RSS2 = (np.array(newTable['y'][2 * nc:]) - X.T.dot(Q_T).T)
  116.     RSS2 = float(RSS2.dot(RSS2.T))
  117.     value = f.ppf(0.95, (n - nc - 2 * k) / 2, (n - nc - 2 * k) / 2)
  118.     if RSS2 / RSS1 < value:
  119.         print('Гипотеза о гомоскедастичности принимается.')
  120.     else:
  121.         print('Гипотеза о гомоскедастичности отвергается.')
  122.  
  123. # 3. Провести оценивание параметров регрессионной модели по доступному обобщенному МНК и по обыкновенному МНК.
  124. # ОМНК
  125. def mnk(u, n, m, p, x1, x2):
  126.     y, x, theta_est, sigma_est = simulation(u, n, m, p, x1, x2)
  127.     X_new = []
  128.     for i in range(len(theta_est)):
  129.         X_new.append([])
  130.         for j in range(n):
  131.             X_new[i].append(f(x1[j], x2[j])[i])
  132.     X_new = np.matrix(X_new).T
  133.     V = np.matrix(np.diag(known_var(y,n)))
  134.     Q_1 = (X_new.T).dot(np.linalg.inv(V))
  135.     Q_1 = np.linalg.inv(Q_1.dot(X_new))
  136.     Q_1 = Q_1.dot((X_new.T).dot(np.linalg.inv(V)))
  137.     Q_1 = Q_1.dot(np.matrix(y).T)
  138.     Q_1 = np.array(Q_1)
  139.     Q_=1
  140.     print('Истинные значения', theta_est)
  141.     #print('МНК', Q_.T) нужно вывести из функции выше
  142.     print('ОМНК', Q_1.T)
  143.     ## 4. Сравнить эффективность оценок в этих двух случаях по квадрату их расстояния до известных истинных значений параметров.
  144.     print('Сравнение МНК с истиной', (theta_est - Q_.T).dot((theta_est - Q_.T).T))
  145.     print('Сравнение ОМНК с истиной', (theta_est - Q_1.T).dot((theta_est - Q_1.T).T))
  146.  
  147.  
  148. def main():
  149.     a = 1
  150.     b = 3
  151.     p = 0.05
  152.     n = 120
  153.     m = 5
  154.     theta = np.array([0.001, 1, 0.001, 0.0001, 1])
  155.     x1, x2 = point_generation(a, b, n)
  156.     u = func_u(x1, x2, theta)
  157.     heteroscedaticity(u, n)
  158.     Breusha_Pagana(u, n)
  159.     #Goldfeld_Quandton(u, n, m, p, x1, x2)
  160.     #mnk(u, n, m, p, x1, x2)
  161.  
  162. if __name__ == '__main__':
  163.     main()
  164.  
RAW Paste Data