<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Least Square Estimation of AR Models and Whitening
Consider an $p$-order autoregressive model driven by white noise.
$$ y_n = \sum_{k=1}^{p} a_ky_{n-k} + \epsilon_n $$
where, $\epsilon$ is a Gaussian white noise with zero mean and some variance $\sigma_{\epsilon}^2$.
Let us assume that we have $N$ sample of this AR process and we are interested in estimating the parameter $a_k$. We can arrange the data into a set of $N-p$ linear equations,
$$ \begin{bmatrix}
y_{0} & y_{1} & \cdots & y_{p-1} \
y_{1} & y_{2} & \cdots & y_{p} \
y_{2} & y_{3} & \cdots & y_{p+1} \
% y_{p+1} & y_{p} & \cdots & y_{2} \
\vdots & \vdots & \ddots & \vdots \
y_{N-p-1} & y_{N-p} & \cdots & y_{N-2}
\end{bmatrix} \begin{bmatrix}
a_p \
a_{p-1} \
\vdots \
a_1
\end{bmatrix} = \begin{bmatrix}
y_{p}\
y_{p+1}\
y_{p+2}\
\vdots\
y_{N-1}
\end{bmatrix} \longrightarrow \mathbf{Y}\mathbf{a} = \mathbf{y} $$
The least-squares solution to this problem is the following,
$$ \hat{\mathbf{a}} = \left(\mathbf{Y}^T\mathbf{Y}\right)^{-1}\mathbf{Y}^T\mathbf{y} $$
Post-multiplying $\mathbf{y}$ by the pseudo-inverse will provide the least square estimate of $\mathbf{a}$.
Whitening using $\hat{a}$
Once $\hat{a}$ is obtained, then the signal $y_n$ can be whitened by passing it through the following moving average filter.
$$w_n = y_n - \sum_{k=1}^{p} \hat{a}ky{n-k}$$
$w_n$ would be the out of this moving average filter, and $w_n$ will be a white noise. The input to this moving average filter the measured time series $y_n$.
Define a stationary AR process & generate a time-series sample
Step1:  Estimate AR parameters using the entire dataset
Step2:  Whiten the AR process output using the  $\hat{a}$
Step3:  Non-stationary AR process with input with varying variance
Even if the input signal was white but non-stationary, the input could be recovered through the AR estimation process, as demonstrated in the following example.
Step4:  Non-stational Autoregressive Model with Varying Model Parameters
Consider the following time-varying $p$-order autoregressive model driven by white noise.
$$ y_n = \sum_{k=1}^{p} a_{n, k}y_{n-k} + \epsilon_n $$
where, $\epsilon$ is a Gaussian white noise with zero mean and some variance $\sigma_{\epsilon}^2$.
Assuming that the parameter variations occur at a larger timescale than the sampling interval, we can estimate these time-varying parameters by carrying out the estimation procedure on a running window of data.
At given point in time $n$, we use the data from the last $M$ samples, ${y_k}{k=n-M+1}^n$. These $M$ data points are used to estimate the AR parameters through the same procedure that was described earlier,
$$ \begin{bmatrix}
y{n-M+1-p} & y_{n-M+2-p} & \cdots & y_{n-M} \
% y_{n-M} & y_{n-M-1} & \cdots & y_{n-M+1-p} \
y_{n-M+2-p} & y_{n-M+3-p} & \cdots & y_{n-M+1} \
% y_{n-M+1} & y_{n-M} & \cdots & y_{n-M+2-p} \
y_{n-M+3-p} & y_{n-M+4-p} & \cdots & y_{n-M+2} \
% y_{n-M+2} & y_{n-M+1} & \cdots & y_{n-M+3-p} \
\vdots & \vdots & \ddots & \vdots \
y_{n-p} & y_{n-p+1} & \cdots & y_{n-1}
\end{bmatrix} \begin{bmatrix}
a_{n,p} \
a_{n,p-1} \
\vdots \
a_{n,1}
\end{bmatrix} = \begin{bmatrix}
y_{n-M+1}\
y_{n-M+2}\
y_{n-M+3}\
\vdots\
y_{n}
\end{bmatrix}$$
$$
\mathbf{Y}{n}\mathbf{a}_n = \mathbf{y}{n,M}
$$
where, $\mathbf{y}{n,M}$ is a column vector consisting of the past $M$ data points, $y_n$ as the last element of the $\mathbf{y}{n,M}$.
$$ \hat{\mathbf{a}}{n} = \left(\mathbf{Y}_n^T\mathbf{Y}_n\right)^{-1}\mathbf{Y}_n^T\mathbf{y}{n, M} $$
Whitening using $\hat{a}_n$
Once $\hat{a}n$ is obtained, then the signal $y_n$ can be whitened by passing it through the following moving average filter.
$$w_n = y_n - \sum{k=1}^{p} \hat{a}{n, k}y{n-k}$$
$w_n$ would be the out of this moving average filter, and $w_n$ will be a white noise. The input to this moving average filter the measured time series $y_n$.
Recursive estimation of AR parameters.
Instead of having to solve a new least squares problem at each time instant, whih can be computaionally intensive, we could estimate $\hat{\mathbf{a}}n$ using $\hat{\mathbf{a}}{n-1}$ and the new data point $y_n$.
The matrix $\mathbf{Y}n$ can be written in terms as the following,
$$ \mathbf{Y}_n = \begin{bmatrix}
\mathbf{y}{n-M, p}^T \
\mathbf{y}{n-M+1, p}^T \
\vdots \
\mathbf{y}{n-1, p}^T
\end{bmatrix} \implies \mathbf{Y}n^T = \begin{bmatrix}
\mathbf{y}{n-M, p} & \mathbf{y}{n-M+1, p} &
\ldots & \mathbf{y}{n-1, p}
\end{bmatrix} $$
$$ \mathbf{Y}n^T\mathbf{Y}_n = \sum{i=1}^{M} \mathbf{y}{n-i, p}\mathbf{y}{n-i, p}^T$$
We can write a recursive expression for $\mathbf{Y}_n^T\mathbf{Y}_n$ as the following,
$$ \mathbf{Y}n^T\mathbf{Y}_n = \mathbf{Y}{n-1}^T\mathbf{Y}{n-1} + \mathbf{y}{n-1, p}\mathbf{y}{n-1, p}^T - \mathbf{y}{n-M-1, p}\mathbf{y}_{n-M-1, p}^T $$
Now, the right hand side can be written as the following,
$$ \mathbf{Y}{n}^T\mathbf{y}{n} = \begin{bmatrix} 
\mathbf{y}{n-M, p} & \mathbf{y}{n-M+1, p} &
\ldots & \mathbf{y}{n-1, p}
\end{bmatrix} \begin{bmatrix}
y{n-M+1}\
y_{n-M+2}\
y_{n-M+3}\
\vdots\
y_{n}
\end{bmatrix} = \sum_{i=1}^{M} y_{n-i+1}\mathbf{y}_{n-i, p}$$
Since, $ \mathbf{Y}{n-1}^T\mathbf{y}{n-1} = \sum_{i=1}^{M} y_{n-i}\mathbf{y}_{n-i-1, p}$, we have
$$ \mathbf{Y}{n}^T\mathbf{y}{n} = \mathbf{Y}{n-1}^T\mathbf{y}{n-1} + y_{n}\mathbf{y}{n-M+1} - y{n-M}\mathbf{y}_{n-M-1, p}$$
We can avoid calculating $\mathbf{Y}n^T\mathbf{Y}_n$ and $\mathbf{Y}_n^T\mathbf{y}{n, M}$ from scratch by using the recursive relationships. Once we calculate these, we can estimate the AR parameters at time  $n$,
$$ \hat{\mathbf{a}}{n} = \left(\mathbf{Y}_n^T\mathbf{Y}_n\right)^{-1}\mathbf{Y}_n^T\mathbf{y}{n, M} $$
<ASSISTANT_TASK:>
Python Code:
# AR parameters
p = 4
a = 1.0 * np.random.rand(p) - 0.5
print "Original AR parameters:\n", a
# Time series data
N = 1000
n = np.arange(0, N)
# Input white noise
eparam = (0, 1.0)
e = np.sqrt(eparam[1]) * np.random.randn(N) + eparam[0]
# Generate AR time series.
y = genARProcess(p, a, eparam, e)
Explanation: Least Square Estimation of AR Models and Whitening
Consider an $p$-order autoregressive model driven by white noise.
$$ y_n = \sum_{k=1}^{p} a_ky_{n-k} + \epsilon_n $$
where, $\epsilon$ is a Gaussian white noise with zero mean and some variance $\sigma_{\epsilon}^2$.
Let us assume that we have $N$ sample of this AR process and we are interested in estimating the parameter $a_k$. We can arrange the data into a set of $N-p$ linear equations,
$$ \begin{bmatrix}
y_{0} & y_{1} & \cdots & y_{p-1} \
y_{1} & y_{2} & \cdots & y_{p} \
y_{2} & y_{3} & \cdots & y_{p+1} \
% y_{p+1} & y_{p} & \cdots & y_{2} \
\vdots & \vdots & \ddots & \vdots \
y_{N-p-1} & y_{N-p} & \cdots & y_{N-2}
\end{bmatrix} \begin{bmatrix}
a_p \
a_{p-1} \
\vdots \
a_1
\end{bmatrix} = \begin{bmatrix}
y_{p}\
y_{p+1}\
y_{p+2}\
\vdots\
y_{N-1}
\end{bmatrix} \longrightarrow \mathbf{Y}\mathbf{a} = \mathbf{y} $$
The least-squares solution to this problem is the following,
$$ \hat{\mathbf{a}} = \left(\mathbf{Y}^T\mathbf{Y}\right)^{-1}\mathbf{Y}^T\mathbf{y} $$
Post-multiplying $\mathbf{y}$ by the pseudo-inverse will provide the least square estimate of $\mathbf{a}$.
Whitening using $\hat{a}$
Once $\hat{a}$ is obtained, then the signal $y_n$ can be whitened by passing it through the following moving average filter.
$$w_n = y_n - \sum_{k=1}^{p} \hat{a}ky{n-k}$$
$w_n$ would be the out of this moving average filter, and $w_n$ will be a white noise. The input to this moving average filter the measured time series $y_n$.
Define a stationary AR process & generate a time-series sample
End of explanation
# Estimate autogressive model parameters
ahat = estARParam(p, y)
# Generate time series using estimated parameters
yhat = genARProcess(p, ahat, eparam, e)
plot_arfit(n, e, y, yhat)
print "Original AR parameters:\n", a
print "Estimated AR parameters:\n", ahat
Explanation: Estimate AR parameters using the entire dataset
End of explanation
w = whiten(y, ahat)
figure(figsize=(3, 3))
plot(e, w, '.')
xlabel("$e_n$", fontsize=16)
ylabel("$w_n$", fontsize=16)
title("$e_n$ vs. $w_n$", fontsize=16)
xticks(fontsize=16)
yticks(fontsize=16);
Explanation: Whiten the AR process output using the  $\hat{a}$
End of explanation
# Input white noise
N = 1000
n = np.arange(0, N)
# v acts as a window that changes the variance of e
v = 1 + (n >= 300) * (n < 700)
e = v * np.sqrt(eparam[1]) * np.random.randn(N) + eparam[0]
# Generate AR time series.
y = genARProcess(p, a, eparam, e)
ahat = estARParam(p, y)
yhat = genARProcess(p, ahat, eparam, e)
plot_arfit(n, e, y, yhat)
print "Original AR parameters:\n", a
print "Estimated AR parameters:\n", ahat
w = whiten(y, ahat)
figure(figsize=(3, 3))
plot(e, w, '.')
xlabel("$e_n$", fontsize=16)
ylabel("$w_n$", fontsize=16)
title("$e_n$ vs. $w_n$", fontsize=16)
xticks(fontsize=16)
yticks(fontsize=16);
Explanation: Non-stationary AR process with input with varying variance
Even if the input signal was white but non-stationary, the input could be recovered through the AR estimation process, as demonstrated in the following example.
End of explanation
%pylab inline
import seaborn as sns
# Stationary AR process implementation
def genARProcess(p, a, eparam, e):
    N = len(e)
    me, ve = eparam
    # Generate the filtered signal
    y = np.zeros(N)
    for i in xrange(p, N):
        y[i] = e[i]
        for k in xrange(p):
            y[i] += a[k] * y[i - k - 1]
    return y
# Estimation of AR parameters for a stationary process,
# using the entire data set.
def estARParam(p, y):
    N = len(y)
    ymat = np.zeros((N - p, p))
    yb = np.zeros((N - p - 1, 1))
    for c in xrange(p, 0, -1):
        ymat[:, p-c] = y[p-c:-c]
    yb = y[p:]
    return np.matmul(np.linalg.pinv(ymat), yb)[::-1]
# Implementation of stationary whitening filter.
def whiten(y, ahat):
    N = len(y)
    p = len(ahat)
    w = np.zeros(N)
    for i in xrange(N):
        w[i] = y[i]
        for k in xrange(p):
            if i - k - 1 >= 0:
                w[i] += - ahat[k] * y[i - k - 1]
    return w
def plot_arfit(n, e, y, yhat):
    figure(figsize=(12, 6))
    subplot(221)
    plot(n, e, label="$e_n$")
    title("$e_n$", fontsize=25)
    xticks(fontsize=18)
    yticks(fontsize=18)
    subplot(223)
    plot(0,0)
    plot(n, y, label="$y_n$")
    xlabel("$n$", fontsize=25)
    title("$y_n$", fontsize=25)
    xticks(fontsize=18)
    yticks(fontsize=18)
    subplot(122)
    plot(y, e, '.', label="$e_n$")
    plot(y, yhat, '.', label="$\hat{y}_n$")
    legend(loc=2, fontsize=25)
    xlabel("$y_n$", fontsize=25)
    ylabel("$e_n, \hat{y}_n$", fontsize=25)
    title("$y_n$ vs. $e_n, \hat{y}_n$", fontsize=25)
    xticks(fontsize=18)
    yticks(fontsize=18)
    tight_layout();
from IPython.core.display import HTML
def css_styling():
    styles = open("styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
Explanation: Non-stational Autoregressive Model with Varying Model Parameters
Consider the following time-varying $p$-order autoregressive model driven by white noise.
$$ y_n = \sum_{k=1}^{p} a_{n, k}y_{n-k} + \epsilon_n $$
where, $\epsilon$ is a Gaussian white noise with zero mean and some variance $\sigma_{\epsilon}^2$.
Assuming that the parameter variations occur at a larger timescale than the sampling interval, we can estimate these time-varying parameters by carrying out the estimation procedure on a running window of data.
At given point in time $n$, we use the data from the last $M$ samples, ${y_k}{k=n-M+1}^n$. These $M$ data points are used to estimate the AR parameters through the same procedure that was described earlier,
$$ \begin{bmatrix}
y{n-M+1-p} & y_{n-M+2-p} & \cdots & y_{n-M} \
% y_{n-M} & y_{n-M-1} & \cdots & y_{n-M+1-p} \
y_{n-M+2-p} & y_{n-M+3-p} & \cdots & y_{n-M+1} \
% y_{n-M+1} & y_{n-M} & \cdots & y_{n-M+2-p} \
y_{n-M+3-p} & y_{n-M+4-p} & \cdots & y_{n-M+2} \
% y_{n-M+2} & y_{n-M+1} & \cdots & y_{n-M+3-p} \
\vdots & \vdots & \ddots & \vdots \
y_{n-p} & y_{n-p+1} & \cdots & y_{n-1}
\end{bmatrix} \begin{bmatrix}
a_{n,p} \
a_{n,p-1} \
\vdots \
a_{n,1}
\end{bmatrix} = \begin{bmatrix}
y_{n-M+1}\
y_{n-M+2}\
y_{n-M+3}\
\vdots\
y_{n}
\end{bmatrix}$$
$$
\mathbf{Y}{n}\mathbf{a}_n = \mathbf{y}{n,M}
$$
where, $\mathbf{y}{n,M}$ is a column vector consisting of the past $M$ data points, $y_n$ as the last element of the $\mathbf{y}{n,M}$.
$$ \hat{\mathbf{a}}{n} = \left(\mathbf{Y}_n^T\mathbf{Y}_n\right)^{-1}\mathbf{Y}_n^T\mathbf{y}{n, M} $$
Whitening using $\hat{a}_n$
Once $\hat{a}n$ is obtained, then the signal $y_n$ can be whitened by passing it through the following moving average filter.
$$w_n = y_n - \sum{k=1}^{p} \hat{a}{n, k}y{n-k}$$
$w_n$ would be the out of this moving average filter, and $w_n$ will be a white noise. The input to this moving average filter the measured time series $y_n$.
Recursive estimation of AR parameters.
Instead of having to solve a new least squares problem at each time instant, whih can be computaionally intensive, we could estimate $\hat{\mathbf{a}}n$ using $\hat{\mathbf{a}}{n-1}$ and the new data point $y_n$.
The matrix $\mathbf{Y}n$ can be written in terms as the following,
$$ \mathbf{Y}_n = \begin{bmatrix}
\mathbf{y}{n-M, p}^T \
\mathbf{y}{n-M+1, p}^T \
\vdots \
\mathbf{y}{n-1, p}^T
\end{bmatrix} \implies \mathbf{Y}n^T = \begin{bmatrix}
\mathbf{y}{n-M, p} & \mathbf{y}{n-M+1, p} &
\ldots & \mathbf{y}{n-1, p}
\end{bmatrix} $$
$$ \mathbf{Y}n^T\mathbf{Y}_n = \sum{i=1}^{M} \mathbf{y}{n-i, p}\mathbf{y}{n-i, p}^T$$
We can write a recursive expression for $\mathbf{Y}_n^T\mathbf{Y}_n$ as the following,
$$ \mathbf{Y}n^T\mathbf{Y}_n = \mathbf{Y}{n-1}^T\mathbf{Y}{n-1} + \mathbf{y}{n-1, p}\mathbf{y}{n-1, p}^T - \mathbf{y}{n-M-1, p}\mathbf{y}_{n-M-1, p}^T $$
Now, the right hand side can be written as the following,
$$ \mathbf{Y}{n}^T\mathbf{y}{n} = \begin{bmatrix} 
\mathbf{y}{n-M, p} & \mathbf{y}{n-M+1, p} &
\ldots & \mathbf{y}{n-1, p}
\end{bmatrix} \begin{bmatrix}
y{n-M+1}\
y_{n-M+2}\
y_{n-M+3}\
\vdots\
y_{n}
\end{bmatrix} = \sum_{i=1}^{M} y_{n-i+1}\mathbf{y}_{n-i, p}$$
Since, $ \mathbf{Y}{n-1}^T\mathbf{y}{n-1} = \sum_{i=1}^{M} y_{n-i}\mathbf{y}_{n-i-1, p}$, we have
$$ \mathbf{Y}{n}^T\mathbf{y}{n} = \mathbf{Y}{n-1}^T\mathbf{y}{n-1} + y_{n}\mathbf{y}{n-M+1} - y{n-M}\mathbf{y}_{n-M-1, p}$$
We can avoid calculating $\mathbf{Y}n^T\mathbf{Y}_n$ and $\mathbf{Y}_n^T\mathbf{y}{n, M}$ from scratch by using the recursive relationships. Once we calculate these, we can estimate the AR parameters at time  $n$,
$$ \hat{\mathbf{a}}{n} = \left(\mathbf{Y}_n^T\mathbf{Y}_n\right)^{-1}\mathbf{Y}_n^T\mathbf{y}{n, M} $$
End of explanation
<END_TASK>