Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- from statistics import mean
- from scipy.stats import chi2
- import pandas as pd
- from scipy.stats import f
- import matplotlib.pyplot as plt
- def func(x1, x2):
- return 1, 1 / x1, x1, x2 ** 2, x2
- def func_u(x1, x2, theta):
- return theta[0] * func(x1, x2)[0] + theta[1] * func(x1, x2)[1] + theta[2] * func(x1, x2)[2] + theta[3] * \
- func(x1, x2)[3] + \
- theta[4] * func(x1, x2)[4]
- def point_generation(a, b, n):
- rng = np.random.default_rng()
- p1 = (b - a) * rng.random(n) + a
- p2 = (b - a) * rng.random(n) + a
- return p1, p2
- def simulation(u, n, m, p, x1, x2):
- omega = np.dot((u - mean(u)), (u - mean(u)).T) / (n - 1)
- disp = np.math.sqrt(p * omega)
- y = u + np.random.normal(0, disp, n)
- # theta estimate
- x = np.array([np.ones(n), func(x1, x2)[1], func(x1, x2)[2], func(x1, x2)[3], func(x1, x2)[4]])
- theta_est = np.linalg.inv(x.dot(x.T)).dot(x).dot(y)
- # sigma estimate
- u_est = func_u(x1, x2, theta_est)
- e = y - u_est
- sigma_est = np.dot(e, e.T) / (n - m)
- return y, x, theta_est, sigma_est
- # Генерация гетероскедатичности
- def heteroscedaticity(u, n):
- y_gen = []
- for i in range(n):
- y_gen.append(u[i] + np.random.normal(0, np.exp(abs(u[i]))))
- return y_gen
- # Функция для генерации вектора известных переменных
- def known_var(var, n):
- v = []
- for i in range(n):
- v.append(float(np.exp(abs(var[i]))))
- v = np.array(v)
- return v
- # 2. Полученные данные проверить по тестам на наличие гетероскедастичности.
- # 1) Тест Бреуша-Пагана (не работает, вид заебись)
- def Breusha_Pagana(u, n):
- y_t = np.array(heteroscedaticity(u, n))
- # МНК и ОМП
- e_t = y_t - u
- sigma_tilda = sum((e_t ** 2) / n)
- c_t = e_t ** 2 / sigma_tilda
- # Построение регрессии и вычисление ESS
- z = np.array([np.ones(n), known_var(y_t, n)])
- alpha = (np.linalg.inv(z.dot(z.T)).dot(z)).dot(c_t.T)
- c_i = z.T * alpha
- c_mean = c_i.mean()
- ESS = 0
- for i in range(n):
- ESS += (c_i[i] - c_mean) ** 2
- a = ESS / 2
- if a < chi2.ppf(0.95, 2):
- print('Гипотеза о гомоскедастичности принимается.')
- else:
- print('Гипотеза о гомоскедастичности отвергается.')
- # 2) Тест Голдфельда-Квандтона (работает, нужно привести код в человеческий вид)
- def Goldfeld_Quandton(u, n, m, p, x1, x2):
- k = m
- y, x, theta_est, sigma_est = simulation(u, n, m, p, x1, x2)
- s = known_var(y, n)
- values = pd.DataFrame({'x1': x1, 'x2': x2, 'y': y, 'u': u, 'exp(|y|)': s})
- # Сортировка по возрастанию
- values = values.sort_values('exp(|y|)')
- values = values.reset_index(drop=True)
- pd.set_option("display.max_rows", None, "display.max_columns", None)
- nc = int(n / 3)
- # Оценка по первым (n-nc)/2 наблюдений
- newTable = values
- X = []
- for i in range(len(theta_est)):
- X.append([])
- for j in range(nc):
- X[i].append(func(newTable['x1'][j], newTable['x2'][j])[i])
- X = np.matrix(X)
- Q_T = np.linalg.inv(X.dot(X.T))
- Q_T = Q_T.dot(X)
- Q_T = Q_T.dot(np.matrix(newTable['y'][0:nc]).T)
- RSS1 = (np.array(newTable['y'][0:nc]) - X.T.dot(Q_T).T)
- RSS1 = float(RSS1.dot(RSS1.T))
- # Оценка по последним (n-nc)/2 наблюдений
- X = []
- for i in range(len(theta_est)):
- X.append([])
- for j in range(nc):
- X[i].append(func(newTable['x1'][j + 2 * nc], newTable['x2'][j + 2 * nc])[i])
- X = np.matrix(X)
- Q_T = np.linalg.inv(X.dot(X.T))
- Q_T = Q_T.dot(X)
- newTable = newTable.reset_index(drop=True)
- Q_T = Q_T.dot(np.matrix(newTable.y[2 * nc:]).T)
- RSS2 = (np.array(newTable['y'][2 * nc:]) - X.T.dot(Q_T).T)
- RSS2 = float(RSS2.dot(RSS2.T))
- value = f.ppf(0.95, (n - nc - 2 * k) / 2, (n - nc - 2 * k) / 2)
- if RSS2 / RSS1 < value:
- print('Гипотеза о гомоскедастичности принимается.')
- else:
- print('Гипотеза о гомоскедастичности отвергается.')
- # 3. Провести оценивание параметров регрессионной модели по доступному обобщенному МНК и по обыкновенному МНК.
- # ОМНК
- def mnk(u, n, m, p, x1, x2):
- y, x, theta_est, sigma_est = simulation(u, n, m, p, x1, x2)
- X_new = []
- for i in range(len(theta_est)):
- X_new.append([])
- for j in range(n):
- X_new[i].append(f(x1[j], x2[j])[i])
- X_new = np.matrix(X_new).T
- V = np.matrix(np.diag(known_var(y,n)))
- Q_1 = (X_new.T).dot(np.linalg.inv(V))
- Q_1 = np.linalg.inv(Q_1.dot(X_new))
- Q_1 = Q_1.dot((X_new.T).dot(np.linalg.inv(V)))
- Q_1 = Q_1.dot(np.matrix(y).T)
- Q_1 = np.array(Q_1)
- Q_=1
- print('Истинные значения', theta_est)
- #print('МНК', Q_.T) нужно вывести из функции выше
- print('ОМНК', Q_1.T)
- ## 4. Сравнить эффективность оценок в этих двух случаях по квадрату их расстояния до известных истинных значений параметров.
- print('Сравнение МНК с истиной', (theta_est - Q_.T).dot((theta_est - Q_.T).T))
- print('Сравнение ОМНК с истиной', (theta_est - Q_1.T).dot((theta_est - Q_1.T).T))
- def main():
- a = 1
- b = 3
- p = 0.05
- n = 120
- m = 5
- theta = np.array([0.001, 1, 0.001, 0.0001, 1])
- x1, x2 = point_generation(a, b, n)
- u = func_u(x1, x2, theta)
- heteroscedaticity(u, n)
- Breusha_Pagana(u, n)
- #Goldfeld_Quandton(u, n, m, p, x1, x2)
- #mnk(u, n, m, p, x1, x2)
- if __name__ == '__main__':
- main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement