%From: Marcilia Andrade Campos \documentstyle[IEEEtran]{article} \begin{document} \pagestyle{myheadings} \tolerance 10000 \markright{APIC'95, El Paso, Extended Abstracts, A Supplement to the international journal of {\rm Reliable Computing}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ } \newcommand{\nat}{\mbox{I${\!}$N}} \newcommand{\real}{\mbox{I${\!}$R}} \title{Interval and High Accuracy Results for Parameter Estimation in Time Series Models} \author{M. A. Campos and M. de B. Correia } \maketitle \begin{abstract} The goal of this paper is to apply interval and high accuracy methods \cite{Moore 1966,Moore 1979,Kulisch and Miranker 1981} to estimating parameters for the AR(1), AR(2), and ARMA(1,1) models \cite{Box and Jenkins 1976}. The proposed interval algorithms are extensions of the Modified Moments Method (MMM) \cite{Bora-Senta and Kounias 1980} to interval values. The results of applying interval methods are compared with the results of applying traditional (floating-point) computations. \end{abstract} \auffil{Departamento de Inform\'atica, Universidade Federal de Pernambuco, Brazil, e-mail:\{mac,mbc\}@di.ufpe.br. The authors are grateful to the referee for his comments.} \section{Introduction} In this paper, by a {\it statistical estimation problem} (or simply an {\it estimation problem}, for short) we mean the following problem: \begin{itemize} \item {\it We know:} \begin{itemize} \item[$\bullet$]a statistical model that describes a certain real-life phenomenon; \item[$\bullet$]the results of experiments and/or observations of this phenomenon. \end{itemize} \item {\it We do not know} the exact values of the parameters that characterize this particular phenomenon. \item {\it We want to estimate} the (unknown) parameters of the model (these parameters are also called {\it population parameters}). \end{itemize} Two kinds of errors are associated with the estimation problem. \begin{itemize} \item The first type of error is a {\it statistical} error. The results of observations form a {\it random sample} from the total population. So, even if all the observations are absolutely precise, instead of the actual value $\mu_p$ of the population parameter, we get an estimate $\mu_s$ that is computed from the random sample values. For example, an estimator for the population mean $ \mu_{p}$ is $ \mu_{s}=(\sum_{i=1}^{n} X_{i})/n$, where $(X_{1},X_{2},\ldots,X_{n})$ is the random sample. \item The second type or error is {\it numerical}. This error is caused by the following two sources: \begin{itemize} \item[$\bullet$] First, the observations and measurements cannot be absolutely accurate. As a result, the observed values $x_i$ are (slightly) different from the actual values $X_i$ of the random sample. \item[$\bullet$]Second, the computations are not absolutely precise. \end{itemize} \item[] After we find an estimator $\mu_{s}$ for the unknown parameter value $\mu_{p}$, we will apply this estimator not to the actual values $X_i$ from the random sample, but to the observed values $x_i$. The resulting estimate $\mu_e$ can thus be different from $\mu_s$. For example, if $(x_{1},x_{2},\ldots,x_{n})$ are the observed values of the random sample, then the estimate for the mean $\mu_{p}$ is $\mu_{e}=(\sum_{i=1}^{n} x_{i})/n$. In general, the estimate $\mu_{s}$ is a function of the random sample, while the estimate $\mu_{e}$ is an {\it evaluation} of the value of that function $\mu_s$ at a given point. \end{itemize} In short: \begin{itemize} \item Finding the function $\mu_s$ that estimates the desired parameter $\mu_p$ is a {\it statistical} problem. \item Estimating the value of this function $\mu_s$ is a {\it numerical} problem. \end{itemize} So, in parameters estimation, there are random and numerical errors. We want to perform our computations in such a way that we will be able to guarantee the precision of the resulting estimate $\mu_e$. The fact that all operations on a computer have a finite precision means that after each operation, we do not have the precise value; if we know the precision $\Delta$ and the computation result $\bar x$, then we can guarantee that the actual value of the result belongs to the interval $[\bar x-\Delta,\bar x+\Delta]$. So, to get operations with guaranteed accuracy, we must perform operations on intervals. A number of methods to define operations on intervals that produce guaranteed precision (in particular, least significant bit accuracy) have been proposed by Moore \cite{Moore 1966,Moore 1979}, Kulisch and Miranker \cite{Kulisch and Miranker 1981}, Alefeld and Herzberger \cite{Alefeld and Herzberger 1983} among others. Implementations of these methods have been developed in extensions of the languages Pascal, Fortran, C and C++ \cite{Adams and Kulisch 1993}. The importance of interval methods for scientific computation is emphasized, for example, by B\"{o}hm in \cite{Bohm 1983}, where the values of the polynomial $f(x)=23616x^{5}-161522x^{4}+401773x^{3}-406754x^{2}+87511x+66576$ are computed for $x=1.7800, 1.7801, ..., 1.7810$. At first, computations are performed using traditional floating-point arithmetic; the results are: $$f(1.7800)=0,\ \ f(1.7801)<0,\ \ f(1.7802)>0,$$ $$f(1.7803)<0, \ \ f(1.7804)>0, \ \ f(1.7805)<0,$$ $$f(1.7806)>0, \ \ f(1.7807)<0, \ \ f(1.7808)<0,$$ $$f(1.7809)>0, \ \ f(1.7810)<0.$$ These results give a misleading information about the location of the real zeroes of this polynomial, because they show that a polynomial of degree 5 has ... 8 zeros. Then, the same values are evaluated by an algorithm based on interval methods and high accuracy arithmetic. The results are: $$f(1.7800)<0, \ \ f(1.7801)<0, \ \ f(1.7802)<0,$$ $$f(1.7803)<0, \ \ f(1.7804)<0, \ \ f(1.7805)>0,$$ $$f(1.7806)>0, \ \ f(1.7807)>0, \ \ f(1.7808)<0,$$ $$f(1.7809)<0, \ \ f(1.7810)<0.$$ The goal of this paper is to apply interval and high accuracy methods to estimating the parameters of the AR(1), AR(2), and ARMA(1,1) time series models. The proposed interval algorithms are extensions of the Modified Moments Method \cite{Bora-Senta and Kounias 1980} to interval values. In the conclusion, we present a comparison between interval and floating-point results. %***** \begin{table*} [htb] \begin{center} \begin{tabular}{ccc|crr} \hline\hline \multicolumn{1}{c}{Series} & \multicolumn{1}{c}{N} & \multicolumn{1}{c}{Model} & \multicolumn{1}{c}{Parameters} & \multicolumn{2}{c}{Parameters Value} \\ \hline\hline & & & & {MMMI} & { MLFP} \\ \cline{5-6} A & 197 & ARMA(1,1) & $\phi_{1}$ & $ 0.8818451_{2}^{3}$ & 0.92\\ & & & & &($\pm$.04) \\ & & & $\theta_{1}$ & $ 0.5143846_{1}^{2}$ & 0.58\\ & & & & &($\pm$.08) \\ \hline C & 226 & AR(1) & $\phi_{1}$ & $ 0.8119850_{7}^{8}$ & 0.82\\ & & & & & ($\pm$.04) \\ \hline D & 310 & AR(1) & $\phi_{1}$ & $ 0.8691036_{7}^{8}$ & 0.87\\ & & & & & ($\pm$.03)\\ \hline E & 100 & AR(2) & $\phi_{1}$ & $ 1.392147_{7}^{8} $ & 1.42\\ & & & & & ($\pm$.07)\\ & & & $\phi_{2}$ & $-0.6945869_{9}^{8}$ & -0.73 \\ & & & & & ($\pm$.07) \\ \hline F & 70 & AR(2) & $\phi_{1}$ & $-0.3065561_{1}^{0}$ & -0.34 \\ & & & & & ($\pm.12$)\\ & & & $\phi_{2}$ & $ 0.1982084_{2}^{3}$ & 0.19 \\ & & & & & ($\pm.12$) \\ \hline\hline \end{tabular} \end{center} \caption{Parameters value of the models AR(1), AR(2) and ARMA(1,1)} \end{table*} %***** Box and Jenkins \cite{Box and Jenkins 1976} define the {\it autoregressive} model of order $p$ (denoted AR($p$)) and the mixed {\it autoregressive-moving average model} of orders $p$ and $q$ (denoted ARMA($p$,$q$)), respectively, by the formulas $$ \phi(B)z_{t}=a_{t},$$ $$ \phi (B)z_{t}=\theta(B)a_{t},$$ where: \begin{itemize} \item $\{ z_{t} \}$, $t\in T=\{0, \pm 1, \pm 2, \ldots \}$ is a sequence that this model describes; \item $\phi(B)=1-\phi_{1}B-\ldots-\phi_{p}B^{p}$ is the {\it autoregressive stationary operator}; \item $\theta(B)=1-\theta_{1}B-\ldots-\theta_{q}B^{q}$ is the {\it moving average invertible operator}; \item $\theta_i$ and $\phi_i$ are real numbers; \item $B$ is the {\it backward shift operator} defined by $B^{m}a_{t}=a_{t-m}$, and \item $\{ a_{t} \}$, $t\in T=\{0, \pm 1, \pm 2, \ldots \}$ is a sequence of uncorrelated random variables with mean zero and constant variance $\sigma^{2}_{a}$. \end{itemize} The modified moments method for autoregressive models was proposed by \cite{Bora-Senta and Kounias 1980}. Its algorithm is based on expressions that use both autocovariance generating function and autocorrelation function (consequently, this method has different versions for AR($p$) and ARMA($p$,$q$) models). In this paper, we propose an interval version of this method: {\it M}odified {\it M}oments {\it M}ethod based on {\it I}ntervals, MMMI for short. \smallskip \noindent{\bf Denotations.} In the following text: \begin{itemize} \item $N$ will denote the number of observed values $z_t$; \item $M$ will denote the number of iterations; \item $m=0,\ldots,M$ will denote the number of the current iteration. \end{itemize} \noindent{\it Comment.} The detailed description of the MMM method can be found in \cite{Bora-Senta and Kounias 1980}. For the readers who are familiar with statistical methods and denotations, we can present the basic iterative formula for MMM: $$\hat{\rho}_{k,m}=$$ $$ (A/\gamma_{0})_{m}(1/N)+(N/(N-k))r_{k}(1-(A/\gamma_{0})_{m}(1/N)),$$ where: \begin{itemize} \item $k=1,\ldots,K$ is the number of parameters to estimate; \item $A=\gamma(1)$, i.e., $A$ is the autocovariance generating function for $B=1$, and \item $r_{k}$ is the sample autocorrelation. \end{itemize} \section{Results} We have designed three interval algorithms and implemented them in Pascal-SC \cite{Rump 1983}. These algorithms compute interval estimates for the parameters of the AR(1), AR(2) and ARMA(1,1) models. In these algorithms, we will distinguish between: \begin{itemize} \item The {\it actual} (unknown) values $\mu_p$ of the parameters; these actual values will be denoted by $\phi_i$ and $\theta_i$. \item The estimates $\mu_s$ of these parameters that are based on the finite sample. These estimates will be denoted by $\tilde\phi_i$ and $\tilde\theta_i$. Since we cannot perform all computations with absolute accuracy, we do not know these values either. \item Instead, we will compute {\it intervals} $I\phi_i$ and $I\theta_i$ that contain the (unknown) values of the estimates $\tilde\phi_i$ and $\tilde\theta_i$. These estimates were previously denoted by $\mu_e$. \end{itemize} \noindent Our three algorithms compute the intervals for the three models: \begin{itemize} \item Algorithm 1 computes the interval $I\phi_1$ that contains the (unknown) precise estimate $\tilde \phi_1$ of the parameter $\phi_{1}$ of the AR(1) model \item[]$z_{t}-\phi_{1}z_{t-1}=a_{t}$. \item Algorithm 2 computes the intervals $I\phi_1$ and $I\phi_2$ that contain the (unknown) precise estimates $\tilde\phi_1$ and $\tilde\phi_2$ of the parameters $\phi_{1}$ and $\phi_{2}$ of the AR(2) model $z_{t}-\phi_{1}z_{t-1}-\phi_{2}z_{t-2}=a_{t}$. \item Algorithm 3 computes the intervals $I\phi_1$ and $I\theta_1$ that contain the (unknown) precise estimates $\tilde\phi_1$ and $\tilde\theta_1$ of the parameters $\phi_{1}$ and $\theta_{1}$ of the ARMA(1,1) model $z_{t}-\phi_{1}z_{t-1}=a_{t}-\theta_{1}a_{t-1}$. \end{itemize} All three algorithms are iterative: on each iteration $m$, we compute the intervals $I\phi_{1,m}$, $I\phi_{2,m}$, and/or $I\theta_{1,m}$ that contain the precise estimates of the corresponding parameters. The intervals that we get on the last iteration form the desired output ($I\phi_1$, $I\phi_2$, or $I\theta_1$). All three algorithms use interval operations. The algorithms were run until one of the following conditions was satisfied: \begin{itemize} \item On the current iteration, we already know the values of the estimates with the desired accuracy $\alpha>0$, i.e., if the width of all the intervals computed on a given iteration does not exceed $2\alpha$ (this number $\alpha$ must be given from the very beginning). \item The number of iterations ($m$) exceeded a given limit $M$ (we chose $M=100$). \item The intervals obtained on two consequent iterations coincided. \end{itemize} Table 1 contains the results of applying our method MMMI to the data from Box and Jenkins \cite{Box and Jenkins 1976}. For comparison, the same table also contains the results of applying maximum likelihood method with floating point arithmetic (MLFP) to the same time series (these results also come from \cite{Box and Jenkins 1976}). Estimates from the column MLFP represent the {\it statistical} error component: namely, they represent the estimates obtained by the maximum likelihood method, and the standard deviations of these estimates. By combining the estimates $e$ and their standard deviations $\sigma$, we can compute the {\it confidence intervals} $(e-k\sigma,e+k\sigma)$ within which the (unknown) actual values $\phi_1$, $\phi_2$, and/or $\theta_1$ lie with a given confidence level (i.e., crudely speaking, with a given probability). For example, if we take a confidence level that correspond to ``one sigma'' deviation, then the numbers corresponding to the series A mean that the actual (unknown) values $\phi_{1}$ and $\theta_{1}$ are inside the confidence intervals (0.88,0.96) and (0.50,0.66) (with this confidence level). Estimates from the column MMMI represent the {\it numerical error} component: for example, the numbers corresponding to series A mean that the (unknown) precise statistical estimates $\tilde \phi_{1}$ and $\tilde \theta_{1}$ for the parameters $\phi_{1}$ and $\theta_{1}$ are inside the intervals $[0.88184512, 0.88184513]$ and $[0.51438461, 0.51438462]$. In particular, if we use a midpoint of each interval as the estimate, then, as a fit to series A, we get an AR(1,1) model $z_{t}-0.88184512z_{t-1}=a_{t}-0.51438461a_{t-1}$. \smallskip \section{Conclusions} Numerical and confidence intervals are {\it different}; their definitions are different, so it is important not to confuse them and not to use the numerical intervals instead the confidence ones. Statistical errors are usually much larger than numerical ones, so, in many applications (especially if we are only interested in confidence intervals), we can safely neglect the numerical errors and consider only statistical ones. However, for some applications, we need the use the statistical estimates themselves; in this case, it is necessary to estimate them as accurately as possible. Our experiments has shown if we use $\alpha\le 10^{-k}$ in the stopping rule, then we get a $k$ digit accuracy in the resulting estimate. For $k=11$ and $k=12$, we have noticed that the larger the size $N$ of the sample, the fewer iterations are needed for our algorithms to convergence, and the smaller are the widths of the resulting intervals. As a reasonable direction of future research in this area, we propose to apply interval and other high accuracy algorithms to improve the accuracy with which the (endpoints of the) confidence intervals are computed. %\bibliography{eq} %\bibliographystyle{alpha} \begin{thebibliography}{99} \bibitem{Adams and Kulisch 1993} E. Adams and U. W. Kulisch, {\bf Scientific Computing with Automatic Result Verification}, Number 189 in Mathematics in Science and Engineering, Academic Press, 1993. \bibitem{Alefeld and Herzberger 1983} G. Alefeld and J. Herzberger, {\bf Introduction to Interval Computation}, Academic Press, N. Y., 1983. \bibitem{Bohm 1983} H. B\"{o}hm,``Evaluation of Arithmetic Expressions with Maximum Accuracy'', {\it A New Approach to Scientific Computation}, 1983, pp. 121--137. \bibitem{Bora-Senta and Kounias 1980} E. Bora-Senta and S. Kounias, ``Parameter Estimation and Order Determination of Autoregressive Models'', {\it Analysing Time Series}, 1980, pp. 93--108. \bibitem{Box and Jenkins 1976} G. E. Box and G. M. Jenkins, {\bf Time Series Analysis Forecasting and Control}, Holden-Day, 1976. \bibitem{Kulisch and Miranker 1981} U. W. Kulisch and W. L. Miranker, {\bf Computer Arithmetic in Theory and Practice}, Academic Press, 1981. \bibitem{Moore 1966} R. E. Moore, {\bf Interval Analysis}, Prentice Hall, Englewood Cliffs, NJ, 1966. \bibitem{Moore 1979} R. E. Moore, {\bf Methods and Applications of Interval Analysis}, SIAM, Philadelphia, Pennsylvania, 1979. \bibitem{Rump 1983} S. M. Rump, ``Computer Demonstration Packages for Standard Problems of Numerical Mathematics'', {\it A New Approach to Scientific Computation}, 1983, pp. 28--49. \end{thebibliography} \end{document}