def fread (filename, sep=';') :
    '''
    Чтение исходных данных из файла CSV или TSV.
    
    Аргументы:
    filename - имя файла с исходными данными
    sep - разделитель ячеек
    
    Результат:
    X, Y - исходные числовые данные
    '''
    f = open(filename, 'r')
    head = f.readline()
    X = []
    Y = []
    for line in f :
        line = line.replace(',', '.').split(sep)
        row = [float(cell) for cell in line[1:]]
        Y.append(row[0])
        X.append([1] + row[1:])
    f.close()
    Y = vector(RDF, Y)
    X = matrix(RDF, X)
    return Y, X
#
def MNK (X, Y) :
    '''
    Построение линейной многофакторной модели методом наименьших квадратов.
    
    Аргументы:
    X, Y - исходные данные
    
    Результат:
    beta_hat - вектор коэффициентов ЛММ
    u_hat - вектор остатков ЛММ
    '''
    beta_hat = (X.transpose()*X)^(-1)*(X.transpose()*Y)
    Y_hat = X * beta_hat
    u_hat = Y - Y_hat
    return beta_hat, u_hat

# график остатков в зависимости от времени
Y, X = fread('data.tsv', '\t')
beta, u = MNK(X, Y)
list_plot(u, size=50, axes_labels=['T', 'residuals'])

# график остатков u_t от u_{t-1}
resid = [(u[i-1], u[i]) for i in xrange(1, len(u))]
list_plot(resid, size=50, color='red', axes_labels=['$u_{t-1}$', '$u_t$'])

def DW (X, Y) :
    '''
    Тест Дарбина - Уотсона на автокорреляцию остатков.
    Метод Кохрейна - Оркатта (расчёт "ро").
    
    Аргументы:
    X, Y - исходные данные
    
    Результат:
    d_st - статистика Дарбина - Уотсона
    rho - оценка параметра "ро" методом Кохрейна - Оркатта
    '''
    beta_hat, u_hat = MNK(X, Y)
    T = len(Y)
    RSS = sum([uh^2 for uh in u_hat])
    d_st = sum([(u_hat[i]-u_hat[i-1])^2 for i in xrange(1, T)])/RSS
    rho = sum([u_hat[i]*u_hat[i-1] for i in xrange(1, T)])/RSS
    return d_st, rho

Y, X = fread('data.tsv', '\t')
T = X.nrows()
N = X.ncols()
# средние значения
Y_mean = mean(Y)
X_mean = vector([mean(X.column(j)) for j in xrange(N)])
# запускаем тест Дарбина - Уотсона
d_st, rho = DW(X, Y)
print 'd_st:', d_st
# вычисляем матрицу преобразования
TA = []
TA.append([sqrt(1-rho^2)]+[0]*(T-1))
for i in xrange(T-1) :
    TA.append([0]*i+[-rho,1]+[0]*(T-i-2))
TA = matrix(RDF, TA)
# преобразуем исходные данные
NewX, NewY = TA*X, TA*Y
# строим новую модель
beta_new = MNK(NewX, NewY)[0]
print 'beta_new:', beta_new
# тест Дарбина - Уотсона для новой модели
print 'd_st_new:', DW(NewX, NewY)[0]

# графики квадратов остатков в зависимости от значений каждой объясняющей переменной
u_new = MNK(NewX, NewY)[1]
u_squared = [p^2 for p in u_new]
for k in xrange(1, N) :
    list_plot(zip(NewX.column(k), u_squared), size=50, color='black', axes_labels=['$x_'+str(k+1)+'$', 'residuals'])
    print

# график квадратов остатков в зависимости от значений регрессанда
list_plot(zip(NewY, u_squared), size=50, axes_labels=['$y$', 'residuals'])

def sortM (M, j0) :
    '''
    Сортировка матрицы M по столбцу j0
    '''
    for i in range (M.nrows()-1) :
        for j in range (M.nrows()-1-i) :
            if M[j,j0] > M[j+1,j0] :
                M.swap_rows(j, j+1)
    return M
#
def GQ (X, Y) :
    '''
    Тест Голдфелда - Куандта на гетероскедастичность остатков
    '''
    XY = matrix(RR, len(Y), 1, Y).augment(X)
    T = XY.nrows()
    N = X.ncols()
    m = round(4/15*T)
    k = T1 = T2 = (T-m)//2
    import scipy.stats
    F_crit = scipy.stats.f.ppf(1-0.1, T2-N, T1-N)
    for j in range (2, XY.ncols()) :
        XY = sortM(XY, j)
        Y1 = XY[:k].column(0)
        X1 = XY[:k,1:]
        beta1, u1 = MNK(X1, Y1)
        Y2 = XY[-k:].column(0)
        X2 = XY[-k:,1:]
        beta2, u2 = MNK(X2, Y2)
        F_stat = sum([x^2 for x in u2])/sum([x^2 for x in u1])*(T1-N)/(T2-N)
        print 'j:', j, 'F_st:', F_stat
        if F_stat > F_crit :
            print 'Гетероскедастичность по x_'+str(j)

GQ(NewX, NewY)

# устранение гетероскедастичности
beta, u = MNK(NewX, NewY)
# строим модель (*) с откликом u^2
u2 = vector([x^2 for x in u])
b, eps = MNK(NewX, u2)
u2_model = NewX*b
# проверяем значимость модели (*)
F_stat = u2_model.dot_product(u2_model)/eps.dot_product(eps)*(T-NewX.ncols())/(NewX.ncols()-1)
import scipy.stats
F_crit = scipy.stats.f.ppf(1-0.1, N-1, T-N)
if F_stat > F_crit :
    print 'Модель (*) значима, что подтверждает гетероскедастичность'
else :
    print 'Модель (*) НЕ значима'
# строим матрицу преобразования
TH = matrix(RR, T, T, [0 for x in range(T^2)])
for i in range (T) :
    TH[i,i] = 1/sqrt(abs(u2_model[i]))
# преобразуем данные
XX, YY = TH*NewX, TH*NewY
# строим новую модель
beta, u = MNK(XX, YY)
print 'beta:', beta
# снова запускаем тест Голдфелда - Куандта
GQ(XX, YY)

def FG (X) :
    '''
    Тест Фаррара - Глобера на мультиколлинеарность в массиве Х
    '''
    T = X.nrows()
    N = X.ncols()
    # нормализация данных
    X_norm = matrix([[(X[i,j]-mean(X.column(j)))/std(list(X.column(j))) for j in xrange(N)] for i in xrange(T)])
    # вычисление корреляционной матрицы
    r = (1/(T-1))*X_norm.transpose()*X_norm
    # хи-квадрат критерий
    chi2_stat = -(T-1-(2*N+5)/6)*ln(det(r))
    import scipy.stats
    chi2_crit = scipy.stats.chi2.isf(0.1, N*(N-1)/2)
    if chi2_stat > chi2_crit :
        print 'В массиве X обнаружена мультиколлинеарность'
    else :
        print 'В массиве X НЕ обнаружена мультиколлинеарность'
    # F-критерий
    Z = r^(-1)
    F_crit = scipy.stats.f.isf(0.1, N-1, T-N)
    for k in xrange(N) :
        F_stat = (Z[k,k]-1)*(T-N)/(N-1)
        if F_stat > F_crit :
            print str(k+2)+'-я переменная мультиколлинеарна с массивом остальных'
    # t-критерий
    t_crit = scipy.stats.t.isf(0.05, T-N)
    for k in xrange(N) :
        for i in xrange(k+1, N) :
            r_ki = -Z[k,i]/sqrt(Z[k,k]*Z[i,i])
            t_stat = abs(r_ki) * sqrt((T-N)/(1-r_ki^2))
            if t_stat > t_crit :
                print 'Обнаружена мультиколлинеарность между x_'+str(k+2)+' и x_'+str(i+2)+': r='+str(r_ki.n(digits=5))

X = fread('data.tsv', '\t')[1][:,1:] # убираем столбец из единиц
FG(X)

def PCA (X) :
    '''
    Анализ главных компонент в массиве X.
    
    Аргумент:
    X - массив исходных данных
    
    Результаты:
    Ip - значения критерия информативности при p = 1, 2, ...
    A - матрица нагрузок главных компонент
    '''
    T = X.nrows()
    N = X.ncols()
    # нормализация данных
    X_norm = matrix([[(X[i,j]-mean(X.column(j)))/std(list(X.column(j))) for j in xrange(N)] for i in xrange(T)])
    # вычисление корреляционной матрицы
    r = (1/(T-1))*X_norm.transpose()*X_norm
    # собственные значения и векторы
    eig = r.eigenvectors_right()
    eig_sorted = sorted(eig, reverse=True)
    eig_values = [p[0] for p in eig_sorted]
    # критерий информативности
    Ip = [(sum(eig_values[:k+1])/sum(eig_values)).n(digits=5) for k in xrange(len(eig_values))]
    # матрица нагрузок
    L = matrix([p[1][0] for p in eig_sorted])
    Lambda = diagonal_matrix([sqrt(p) for p in eig_values])
    A = L.transpose()*Lambda
    return Ip, A

Ip, A = PCA(X)
show('$I_p='+latex(Ip)+'$')
show('$A='+latex(A)+'$')