Advertisement
evgeniya_polyntseva

smad 4

Dec 1st, 2020
669
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 6.13 KB | None | 0 0
  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.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement