%From: Marcilia Andrade Campos <mac@di.ufpe.br>

\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}
 

