data = sio.loadmat('Tut3_file1.mat')
DLPFC = data['DLPFC']
x = DLPFC[0, :]
plt.plot(DLPFC.T)

x = DLPFC[0]
X_p, x_T = set_data(3, x)

res = LA.lstsq(X_p, x_T)

eps = x_T - np.dot(X_p, res[0])
XpXp = np.diag(LA.inv(np.dot(X_p.T, X_p)))
sdeps = np.std(eps, axis=0, ddof=0)
res[0].flatten() / (sdeps * np.sqrt(XpXp))
err = (sdeps * np.sqrt(XpXp))

"""
Explanation: Task 1
End of explanation
"""


import statsmodels.api as sm
ols = sm.OLS(x_T, X_p).fit()
ols.summary()

"""
Explanation: The t-statistic computes how distant our parameters are from the zero hypothesis in terms of standard errors. So for $a_0$ and $a_1$ we cannot refute the zero hypothesis, but in case of $a_2$ and $a_3$ it is highly unlikely that  our that comes from random fluctuations.
End of explanation
"""


ols.params, res[0].flatten()
ols.HC1_se, 


"""
Explanation: The likelihood function is:
$ f (Y) = (\frac{1}{\sqrt{2 \pi} \sigma} \exp{- \frac{1}{2 \sigma^2}\sum^n_{i=1}(Y_i - a_0 X_{i1} - ... - a_4)}) $, so the log likelihood will be:
$-n\ln(\sqrt{2 \pi} \sigma) - \frac{1}{2 \sigma^2}\sum^n_{i=1}(Y_i - a_0 X_{i1} - ... - a_4)$.
The result in our case can be seen in the summary of the previous cell.
End of explanation
"""


X = np.ones((5, len(x)))
X[1:3, :] = DLPFC
X[3:, :] = data['Parietal']
X_T = X[1:, 1:].T.reshape((359, 4))
A = np.zeros((4, 5))
tv = np.zeros((4, 5))
for i in range(1, 5):

    X_T = X[i, 1:].T.reshape((359, 1))
    X_p = X[:, :-1].T
    ols = sm.OLS(X_T, X_p).fit()
    A[i-1, :] = ols.params
    tv[i-1, :] = ols.pvalues

"""
Explanation: Task 2
We do the VAR(1) regression similar to the previous case using all the time series.
End of explanation
"""


A[:, :]

X_T = X[1:, 1:].T.reshape((359, 4))

res = LA.lstsq(X_p, X_T)

"""
Explanation: We can see on the off diagonal terms are much smaller than the terms on diagonal. That shows that there is not as much influence between features as on the feature itself.
End of explanation
"""


res[0].shape

"""
Explanation: For checking if it si astationaty process  we will see if the  eignvalues are in the unitcircle. They are.
End of explanation
"""


eps = X_T - np.dot(X_p, res[0])
sdeps = np.std(eps, axis=0)

XpXp = LA.inv(np.dot(X_p.T, X_p))

"""
Explanation: The loglikelihood of this model will be the
End of explanation
"""


