diff --git a/Jenkinsfile b/Jenkinsfile index 43acc44..ac8f2ec 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -12,6 +12,6 @@ node { } stage('Publish PDF') { - sh 'scp -i /root/.ssh/id_rsa ats-zf.pdf thisfro@192.168.178.45:/opt/containers/apache2/html/download/latex-previews/ats-zf.pdf' + sh 'scp -i /root/.ssh/id_rsa ats-zf.pdf jannis@192.168.178.45:/var/www/html/download/latex-previews/ats-zf.pdf' } } \ No newline at end of file diff --git a/img/checkresiduals.png b/img/checkresiduals.png deleted file mode 100644 index 85fcdd5..0000000 Binary files a/img/checkresiduals.png and /dev/null differ diff --git a/img/pacf.png b/img/pacf.png deleted file mode 100644 index 7a5cb64..0000000 Binary files a/img/pacf.png and /dev/null differ diff --git a/main.tex b/main.tex index e1b3ba4..d80367a 100644 --- a/main.tex +++ b/main.tex @@ -180,7 +180,7 @@ mathematically formulated by strict stationarity. \end{tabular} \subsubsection{Weak} \label{weak-stationarity} - It is impossible to «prove» the theoretical concept of stationarity from data. We can only search for evidence in favor or against it. \\ + It is impossible to "prove" the theoretical concept of stationarity from data. We can only search for evidence in favor or against it. \\ \vspace{0.1cm} However, with strict stationarity, even finding evidence only is too difficult. We thus resort to the concept of weak stationarity. @@ -522,346 +522,6 @@ var.ts <- 1/n^2*acf(b,lag=0,type="cov")$acf[1]*(n+2*sum(((n-1):(n-10))*acf(b,10) mean(b) + c(-1.96,1.96)*sqrt(var.ts) \end{lstlisting} -\subsection{Partial autocorrelation (PACF)} -The $k$-th partial autocorrelation $\pi_k$ is defined as the correlation between $X_{t+k}$ and $X_t$, given all the values in between. -$$\pi_k = Cor(X_{t+k},X_t | X_{t+1},...,X_{t+k-1} = x_{t+k-1})$$ -\begin{itemize} - \item Given a time series X t , the partial autocorrelation of lag $k$, is the autocorrelation between $X_t$ and $X_{t+k}$ with the linear dependence of $X_{t+1}$ through to $X_{t+k-1}$ removed. - \item One can draw an analogy to regression. The ACF measures the „simple“ dependence between $X_t$ and $X_{t+k}$, whereas the PACF measures that dependence in a „multiple“ fashion.\footnote{See e.g. \href{https://n.ethz.ch/~jannisp/download/Mathematik-IV-Statistik/zf-statistik.pdf}{\textit{Mathematik IV}}} -\end{itemize} -$$\pi_1 = \rho_1$$ -$$\pi_2 = \frac{\rho_2 - \rho_1^2}{1-\rho_1^2}$$ -for AR(1) moderls, we have $\pi_2 = 0$, because $\rho_2 = \rho_1^2$, i.e. there is no conditional relation between $(X_t, X_{t+2} | X_{t+1})$ - -\begin{lstlisting}[language=R] -pacf(wave, ylim=c(1,1)) -\end{lstlisting} - -\begin{figure}[H] - \centering - \includegraphics[width=.25\textwidth]{pacf.png} - \caption{PACF for wave tank} - \label{fig:pacf} -\end{figure} - -\section{Basics of modelling} -\subsection{White noise} -\begin{quote} - A time series $(W_1, W_2,..., W_n)$ is a \textbf{White Noise} series if the random variables $W_1 , W_2,...$ are i.i.d with mean zero. -\end{quote} -This implies that all $W_t$ have the same variance $\sigma_W^2$ and -$$Cov(W_i,W_j) = 0 \, \forall \, i \neq j$$ -Thus, there is no autocorrelation either: $\rho_k = 0 \, \forall \, k \neq 0$. \\ -\vspace{.2cm} -If in addition, the variables also follow a Gaussian distribution, i.e. $W_t \sim N(0, \sigma_W^2)$, the series is called \textbf{Gaussian White Noise}. The term White Noise is due to the analogy to white light (all wavelengths are equally distributed). - -\subsection{Autoregressive models (AR)} -In an $AR(p)$ process, the random variable $X_t$ depends on an autoregressive linear combination of the preceding $X_{t-1},..., X_{t-p}$, plus a „completely independent“ term called innovation $E_t$. -$$X_t = \alpha_1 X_{t-1} + ... + \alpha_p X_{t-p} + E_t$$ -Here, $p$ is called the order of the AR model. Hence, we abbreviate by $AR(p)$. An alternative notation is with the backshift operator $B$: -$$(1-\alpha_1 B - \alpha_2 B^2 - ... \alpha_p B^p) X_t = E_t \Leftrightarrow \Phi(B)X_t = E_t$$ -Here, $\Phi(B)$ is called the characteristic polynomial of the $AR(p)$. It determines most of the relevant properties of the process. - -\subsubsection{AR(1)-Model}\label{ar-1} -$$X_t = \alpha_1 X_{t-1} + E_t$$ -where $E_t$ is i.i.d. with $E[E_t] = 0$ and $Var(E_t) = \sigma_E^2$. We also require that $E_t$ is independent of $X_s, s p$. The behavior before lag $p$ can be arbitrary. -\end{itemize} -If what we observe is fundamentally different from the above, it is unlikely that the series was generated from an $AR(p)$-process. We thus need other models, maybe more sophisticated ones. - -\subsubsection{Parameter estimation} -Observed time series are rarely centered. Then, it is inappropriate to fit a pure $AR(p)$ process. All R routines by default assume the shifted process $Y_t = m + X_t$. Thus, we face the problem: -$$(Y_t - m) = \alpha_1(Y_{t-1} - m) + ... + \alpha_p(Y_{t-p} - m) + E_t$$ -The goal is to estimate the global mean m , the AR-coefficients $\alpha_1 ,..., \alpha_p$, and some parameters defining the distribution of the innovation $E_t$. We usually assume a Gaussian, hence this is $\sigma_E^2$.\\ -\vspace{.2cm} -We will discuss 4 methods for estimating the parameters:\\ -\vspace{.2cm} - -\textbf{OLS Estimation} \\ -If we rethink the previously stated problem, we recognize a multiple linear regression problem without -intercept on the centered observations. What we do is: -\begin{enumerate} - \item Estimate $\hat{m} = \bar{y}$ and $x_t = y_t - m$ - \item Run a regression without intercept on $x_t$ to obtain $\hat{\alpha_1},\dots,\hat{\alpha_p}$ - \item For $\hat{\sigma_E^2}$, take the residual standard error from the output -\end{enumerate} - -\vspace{.2cm} - -\textbf{Burg's algorithm} \\ -While OLS works, the first $p$ instances are never evaluated as responses. This is cured by Burg’s algorithm, which uses the property of time-reversal in stochastic processes. We thus evaluate the RSS of forward and backward prediction errors: -$$\sum_{t=p+1}^n \bigg[\bigg(X_t - \sum_{k=1}^p \alpha_k X_{t-k}\bigg)^2 + \bigg(X_{t-p} - \sum_{k=1}^p \alpha_k X_{t-p+k}\bigg)^2 \bigg]$$ -In contrast to OLS, there is no explicit solution and numerical optimization is required. This is done with a recursive method called the Durbin-Levison algorithm (implemented in R). - -\begin{lstlisting}[language=R] -f.burg <- ar.burg(llynx, aic=F, order.max=2) -\end{lstlisting} - -\vspace{.2cm} - -\textbf{Yule-Walker Equations} \\ -The Yule-Walker-Equations yield a LES that connects the true ACF with the true AR-model parameters. We plug-in the estimated ACF coefficients: -$$\hat{\rho}(k) = \hat{\alpha_k}\hat{\rho}(k-1) + \dots + \hat{\alpha_p}\hat{\rho}(k-p), \, \mathrm{for} \, k=1,\dots,p$$ -and solve the LES to obtain the AR-parameter estimates.\\ -\vspace{.2cm} -In R we can use \verb|ar.yw()| \\ - -\vspace{.2cm} - -\textbf{Maximum-likelihood-estimation} \\ -Idea: Determine the parameters such that, given the observed time series $(y_1 ,\dots, y_n)$, the resulting model is the most plausible (i.e. the most likely) one. \\ -This requires the choice of a probability model for the time series. By assuming Gaussian innovations, $E_t \sim N (0,\sigma_E^2)$ , any $AR(p)$ process has a multivariate normal distribution: -$$Y = (Y_1,\dots,Y_n) \sim N(m \cdot \vec{1},V)$$ -with $V$ depending on $\vec{\alpha},\sigma_E^2$ \\ -MLE then provides simultaneous estimates by optimizing: -$$L(\alpha,m,\sigma_E^2) \propto \exp \bigg( \sum_{t=1}^n(x_t - \hat{x_t}) \bigg)$$ - -\begin{lstlisting}[language=R] -> f.ar.mle -Call: arima(x = log(lynx), order = c(2, 0, 0)) -\end{lstlisting} - -\vspace{.2cm} - -\textbf{Some remarks} \\ -\begin{itemize} - \item All 4 estimation methods are asymptotically equivalent and even on finite samples, the differences are usually small. - \item All 4 estimation methods are non-robust against outliers and perform best on data that are approximately Gaussian. - \item Function \verb|arima()| provides standard errors for $\hat{m}; \hat{\alpha}_1 ,\dots, \hat{\alpha}_p$ so that statements about significance become feasible and confidence intervals for the parameters can be built. - \item \verb|ar.ols()|, \verb|ar.yw()| and \verb|ar.burg()| allow for convenient choice of the optimal model order $p$ using the AIC criterion. Among these methods, \verb|ar.burg()| is usually preferred. - -\end{itemize} - -\subsection{Model diagnostics} -\subsubsection{Residual analysis}\label{residual-analysis} -"residuals" = "estimated innovations" -$$\hat{E_t} = (y_t - \hat{m}) - (\hat{\alpha_1}(y_{t-1} - \hat{m}) - \dots - \hat{\alpha}_p(y_{t-1} - \hat{m}))$$ -With assumptions as in Chapter \ref{ar-1} \\ - -\vspace{.2cm} -We can check these, using (in R: \verb|tsdisplay(resid(fit))|) -\begin{itemize} - \item Time-series plot of $\hat{E}_t$ - \item ACF/PACF-plot of $\hat{E}_t$ - \item QQ-plot of $\hat{E}_t$ -\end{itemize} - -The time-series should look like white-noise \\ -\vspace{.2cm} -\textbf{Alternative} \\ -Using \verb|checkresiduals()|: \\ -A convenient alternative for residual analysis is this function from \verb|library(forecast)|. It only works correctly when fitting with \verb|arima()|, though. - -\begin{lstlisting}[language=R] -> f.arima <- arima(log(lynx), c(11,0,0)) -> checkresiduals(f.arima) -Ljung-Box test -data: Residuals from ARIMA(11,0,0) with non-zero mean -Q* = 4.7344, df = 3, p-value = 0.1923 -Model df: 12. Total lags used: 15 -\end{lstlisting} - -The function carries out a Ljung-Box test to check whether residuals are still correlated. It also provides a graphical output: -\begin{figure}[H] - \centering - \includegraphics[width=.25\textwidth]{checkresiduals.png} - \caption{Example output from above code} - \label{fig:checkresiduals} -\end{figure} - -\subsubsection{Diagsnostic by simulation} -As a last check before a model is called appropriate, simulating from the estimated coefficients and visually inspecting the resulting series (without any prejudices) to the original one can be beneficial. -\begin{itemize} - \item The simulated series should "look like" the original. If this is not the case, the model failed to capture (some of) the properties in the original data. - \item A larger or more sophisticated model may be necessary in cases where simulation does not recapture the features in the original data. -\end{itemize} - -\subsection{Moving average models (MA)} -Whereas for $AR(p)$-models, the current observation of a series is written as a linear combination of its own past, $MA(q)$-models can be seen as an extension of the "pure" process -$$X_t = E_t$$ -in the sense that the last q innovation terms $E_{t-1} , E_{t-2} ,...$ are included, too. We call this a moving average model: -$$X_t = E_t + \beta_1 E_{t-1} + \beta_2 E_{t-2} + \dots + \beta_q E_{t-q}$$ -This is a time series process that is stationary, but not i.i.d. In many aspects, $MA(q)$ models are complementary to $AR(p)$. - -\subsubsection{Stationarity of MA models} -We first restrict ourselves to the simple $MA(1)$-model: -$$X_t = E_t + \beta_1 E_{t-1}$$ -The series $X_t$ is always weakly stationary, no matter what the choice of the parameter $\beta_1$ is. - -\subsubsection{ACF/PACF of MA processes} -For the ACF -$$\rho(1) = \frac{\gamma(1)}{\gamma(0)} = \frac{\beta_1}{1+\beta_1^2} < 0.5$$ -and -$$\rho(k) = 0 \, \forall \, k > 1$$ - -Thus, we have a «cut-off» situation, i.e. a similar behavior to the one of the PACF in an $AR(1)$ process. This is why and how $AR(1)$ and $MA(1)$ are complementary. - -\subsubsection{Invertibility} -Without additional assumptions, the ACF of an $MA(1)$ does not allow identification of the generating model. -$$X_t = E_t + 0.5 E_{t-1}$$ -$$U_t = E_t + 2 E_{t-1}$$ -have identical ACF! -$$\rho(1) = \frac{\beta_{1}}{1+\beta_1^2} = \frac{1/\beta_1}{1+(1/\beta_1^2)}$$ - -\begin{itemize} - \item An $MA(1)$-, or in general an $MA(q)$-process is said to be invertible if the roots of the characteristic polynomial $\Theta(B)$ exceed one in absolute value. - \item Under this condition, there exists only one $MA(q)$-process for any given ACF. But please note that any $MA(q)$ is stationary, no matter if it is invertible or not. - \item The condition on the characteristic polynomial translates to restrictions on the coefficients. For any MA(1)-model, $|\beta_1| < 1$ is required. - \item R function \verb|polyroot()| can be used for finding the roots. -\end{itemize} - -\textbf{Practical importance:} \\ -The condition of invertibility is not only a technical issue, but has important practical meaning. All invertible $MA(q)$ processes can be expressed in terms of an $AR(\infty)$, e.g. for an $MA(1)$: -\begin{align*} -X_t &= E_t + \beta_1 E_{t-1} \\ - &= E_t + \beta_1(X_{t-1} - \beta_1 E_{t-2}) \\ - &= \dots \\ - &= E_t + \beta_1 X_{t-1} - \beta_1^2 X_{t-2} + \beta_1^3X_{t-3} + \dots \\ - &= E_t + \sum_{i=1}^\infty \psi_i X_{t-i} -\end{align*} - -\subsection{Fitting MA(q)-models to data} -As with AR(p) models, there are three main steps: -\begin{enumerate} - \item Model identification - \begin{itemize} - \item Is the series stationary? - \item Do the properties of ACF/PACF match? - \item Derive order $q$ from the cut-off in the ACF - \end{itemize} - \item Parameter estimation - \begin{itemize} - \item How to determine estimates for $m, \beta_1 ,\dots, \beta_q, \sigma_E^2$? - \item Conditional Sum of Squares or MLE - \end{itemize} - \item Model diagnostics - \begin{itemize} - \item With the same tools/techniques as for AR(p) models - \end{itemize} -\end{enumerate} - -\subsubsection{Parameter estimation} -The simplest idea is to exploit the relation between model parameters and autocorrelation coefficients («Yule-Walker») after the global mean $m$ has been estimated and subtracted. \\ -In contrast to the Yule-Walker method for AR(p) models, this yields an inefficient estimator that generally generates poor results and hence should not be used in practice. - -\vspace{.2cm} -It is better to use \textbf{Conditional sum of squares}:\\ -This is based on the fundamental idea of expressing $\sum E_t^2$ in terms of $X_1 ,..., X_n$ and $\beta_1 ,\dots, \beta_q$, as the innovations themselves are unobservable. This is possible for any invertible $MA(q)$, e.g. the $MA(1)$: -$$E_t = X_t = \beta_1 X_{t-1} + \beta_1^2 X_{t-2} + \dots + (-\beta)^{t-1} X_1 + \beta_1^t E_0$$ -Conditional on the assumption of $E_0 = 0$ , it is possible to rewrite $\sum E_t^2$ for any $MA(1)$ using $X_1 ,\dots, X_n $ and $\beta_1$. \\ -\vspace{.2cm} -Numerical optimization is required for finding the optimal parameter $\beta_1$, but is available in R function \verb|arima()| with: -\begin{lstlisting}[language=R] -> arima(..., order=c(...), method="CSS") -\end{lstlisting} - -\textbf{Maximium-likelihood estimation} -\begin{lstlisting}[language=R] -> arima(..., order=c(...), method="CSS-ML") -\end{lstlisting} -This is the default methods in R, which is based on finding starting values for MLE using the CSS approach. If assuming Gaussian innovations, then: -$$X_t = E_t + \beta_1 E_{t-1} + \beta_q E_{t-q}$$ -will follow a Gaussian distribution as well, and we have: -$$X = (X_1, \dots, X_n) \sim N(0,V)$$ -Hence it is possible to derive the likelihood function and simultaneously estimate the parameters $m;\beta_1,\dots,\beta_q;\sigma_E^2$. - -\subsubsection{Residual analysis} -See \ref{residual-analysis} - -\subsection{ARMA(p,q)-models} -An $ARMA(p,q)$ model combines $AR(p)$ and $MA(q)$: -$$X_t = \alpha_1 X_{t-1} + \dots + \alpha_p X_{t-p} + E_t + \beta_1 E_{t-1} + \dots + \beta_q E{t-q}$$ -where $E_t$ are i.i.d. innovations (=a white noise process).\\ -\vspace{.2cm} -It‘s easier to write $ARMA(p,q)$’s with the characteristic polynomials: \\ -\vspace{.2cm} -$\Phi(B)X_t = \Theta(B)E_t$, where \\ -$\Phi(z) = 1 - \alpha_1 z - \dots - \alpha_p z^p$, is the cP of the $AR$-part, and \\ -$\Theta(z) = 1 + \beta_1 z + \dots + \beta_1 z^q$ is the cP of the $MA$-part - -\subsubsection{Properties of ARMA(p,q)-Models} -The stationarity is determined by the $AR(p)$-part of the model:\\ -If the roots of the characteristic polynomial $\Phi(B)$ exceed one in absolute value, the process is stationary.\\ -\vspace{.2cm} -The invertibility is determined by the $MA(q)$-part of the model:\\ -If the roots of the characteristic polynomial $\Theta(B)$ exceed one in absolute value, the process is invertible.\\ -\vspace{.2cm} -Any stationary and invertible $ARMA(p,q)$ can either be rewritten in the form of a non-parsimonious $AR(\infty)$ or an $MA(\infty)$.\\ -In practice, we mostly consider shifted $ARMA(p,q)$: $Y_t = m + X_t$ - -\begin{table}[H] - \centering - \begin{tabular}{l|l|l} - & ACF & PACF \\ - \hline - $AR(p)$ & exponential decay & cut-off at lag $p$ \\ - $MA(q)$ & cut-off at lag $q$ & exponential decay \\ - $ARMA(p,q)$ & mix decay/cut-off & mix decay/cut-off \\ - \end{tabular} - \caption{Comparison of $AR$-,$MA$-, $ARMA$-models} -\end{table} - -\begin{itemize} - \item In an $ARMA(p,q)$, depending on the coefficients of the model, either the $AR(p)$ or the $MA(q)$ part can dominate the ACF/PACF characteristics. - \item In an $ARMA(p,q)$, depending on the coefficients of the model, either the $AR(p)$ or the $MA(q)$ part can dominate the ACF/PACF characteristics. - -\end{itemize} - -\subsubsection{Fitting ARMA-models to data} -See $AR$- and $MA$-modelling - -\subsubsection{Identification of order (p,q)} -May be more difficult in reality than in theory: -\begin{itemize} - \item We only have one single realization of the time series with finite length. The ACF/PACF plots are not «facts», but are estimates with uncertainty. The superimposed cut-offs may be difficult to identify from the ACF/PACF plots. - \item $ARMA(p,q)$ models are parsimonius, but can usually be replaced by high-order pure $AR(p)$ or $MA(q)$ models. This is not a good idea in practice, however! - \item In many cases, an AIC grid search over all $ARMA(p,q)$ with $p+q < 5$ may help to identify promising models. -\end{itemize} - \scriptsize \section*{Copyright} @@ -872,8 +532,8 @@ Jannis Portmann, FS21 \section*{References} \begin{enumerate} - \item ATSA\_Script\_v210219.docx, M. Dettling - \item ATSA\_Slides\_v210219.pptx, M. Dettling + \item ATSA\_Script\_v219219.docx, M. Dettling + \item ATSA\_Slides\_v219219.pptx, M. Dettling \end{enumerate} \section*{Image sources}