#!/bin/sh
# This is a shell archive (produced by shar 3.49)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#
# made 06/07/1993 14:59 UTC by ikaranta@hutcs.cs.hut.fi
# Source directory /m/huuhaa/users11/ikaranta/tsintser
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length  mode       name
# ------ ---------- ------------------------------------------
#   8511 -rw-r--r-- README.DOC
#  10812 -rw-rw-rw- aep.m
#    908 -rw-r--r-- app.m
#   1672 -rw-r--r-- autocorr.m
#    597 -rw-r--r-- broyden.m
#   1268 -rw-r--r-- forec.m
#   1569 -rw-r--r-- fval.m
#   3429 -rw-r--r-- inter1.m
#   2767 -rw-r--r-- inter1m.m
#   3692 -rw-r--r-- inter1r.m
#   2748 -rw-r--r-- inter2.m
#   4115 -rw-r--r-- inter2r.m
#   1858 -rw-r--r-- intvar.m
#    848 -rw-r--r-- pttopi.m
#
# ============= README.DOC ==============
if test -f 'README.DOC' -a X"$1" != X"-c"; then
	echo 'x - skipping README.DOC (File already exists)'
else
echo 'x - extracting README.DOC (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'README.DOC' &&
Documentation of the time series interpolation package tsinter
--------------------------------------------------------------
X
This document describes the matlab m-functions made to find estimates
for missing observations in a time series. The algorithms are
described in Karanta (1993). It is assumed that an ARIMA model
\begin{equation}
\phi(B)z_t=\theta(B)a_t
\end{equation}
describing the series is known. Here $z_t$ is the value of the time
series variable at time instant $t$, $a_t$ is white noise, $B$ is 
the backward shift operator (that is, $Bz_t=z_{t-1}$). $\phi(B)=
1-\phi_1B-...\phi_pB^p$ and $\theta(B)=1-\theta_1B-...\theta_qB^q$.
X
These programs may be freely distributed and used provided that the
header comments and this file are included. The programs are
distributed as is, with no guarantee on their correctness, functioning
or suitability to any specific purpose.
X
This document covers the following m functions:
X
Primary functions:
inter1.m       interpolation for time series with many missing observations
inter1r.m      same as inter1.m, but contains more user-choosable parameters
inter1m.m      interpolation for time series when information on some parts 
X               of the series is restricted to the mean values of each of 
X               these parts
inter2.m       interpolation for time series with few missing observations
inter2r.m      same as inter2.m, but contains more user-choosable parameters
intvar.m       calculates the error covariance matrix.
forec.m        calculate forecasts for a given lead time with a given
X               time series with no missing observations (actually missing
X               observations have been replaced by estimates) and the
X               noises corresponding to the time series.
X
Secondary functions:
X  app.m        calculate the value of a given time series lag polynomial.
X               Needed by all the interpolation routines.
X  autocorr.m   calculate the autocorrelations of an ARIMA process. Needed by
X               intvar.m
X  broyden.m    calculate the approximate Hessian matrix in Broyden's
X               method. Needed by inter1.m, inter1r.m and inter2.m
X  fval.m       estimates for the time series (actual if known and simulated
X               if unknown) and the noise, and finds the norm of the noise.
X  pttopi.m     calculate the quotient of two lag polynomials. Needed by
X               autocorr.m.
X
For testing purposes, the package includes also the matrix file
aep.m.
X
X
Argument variables of the programs
----------------------------------
X
The order in which the arguments should be fed to the functions 
can be seen from the functions themselves. This information can
also be caught by typing help <function>, e.g., help inter1.
X
The time series z with missing observations should be given as a
real-valued column vector with the value NaN on the place of a
missing observation. Thus, for example, the following vector is
valid:
3.2
2.4 
NaN
5.6
X
It is not recommended that unknown values should exist before the 
first known value or after the last known value.
X
The vector zest needed by the functions inter1r and inter2r should
contain initial estimates of the missing values at their proper
places. Otherwise the values of zest don't matter. It is recommended,
though, that the other values be set to the actual values of the time
series in order to avoid confusion and put the estimates at their
right places.
X
The noise variables a (if required by the program) should be given as
a real-valued column vector. If this is required by the program, it
may well be set to zero if no better information about its values is
available.
X
The phi matrix should contain the \phi parameters together with the
corresponding lags. The zeroth order term 1 is omitted. Hence, the
polynomial $1-0.8B-B+0.1B^{12}$ is represented by the matrix
0.8   1
1     2
-0.1  12
X
The theta matrix should contain the \theta parameters together with the
corresponding lags. Again, the zeroth order term 1 is omitted. The format
is as with the phi matrix.
X
The means matrix needed in the function inter1m should contain the
mean values for the periods of interest together with the
corresponding start and end points of the interval. So, if we are
interested in two intervals for which we have the mean values, say,
the mean of interval 1 being 3.1 and start point 8, end point 14; and
the mean of interval 2 being 4.5 and start point 20, end point 32, we
would have the following means matrix:
3.1    8    14
4.5   20    32
X
The parameter tol (needed in functions inter1r and inter2r) is the
tolerance in the stopping criterion of the algorithms. It is the number
the square of the norm of the gradient vector should be so that the 
iteration can be stopped. In functions not requiring this parameter
as input the tolerance is set to 1E-12.
X
Output variables of the program
-------------------------------
X
zhat is a column vector containing the time series together with 
estimates for the missing observations at appropriate places.
X
ahat is a column vector containing the noise values (calculated by the
functions) corresponding to zhat.
X
Test series: aep.m
--------------------
X
The package includes the matrix file aep.m. It is an ordinary matlab
m-file based on the international airline passenger rates from 1949 to
1960. The series was handled in chapter 9 of Box and Jenkins (1976).
There, the model
X
(1-B)(1-B^{12})z_t=(1-0.396B)(1-0.614B^{12})a_t
X
was obtained for the series.  (Here z_t is the international airline
passenger time series, a_t is white noise, and B is the lag operator).
X
The file contains the following vectors and matrices.
aev   the monthly international airline passenger rates from 1949 to 1960.
X      The series is series J from Box and Jenkins (1976).
ae1   the monthly international airline passenger rates from 1949 to 1960
X      with the observations for January-November for the years 1955-1960
X      missing
ae2   the monthly international airline passenger rates from 1949 to 1960
X      with the observations from July 1949, June 1957, July 1957, 
X      August 1957, and July 1960 missing
ae3   the monthly international airline passenger rates from 1949 to 1960
X      with all July observations missing and June and August 1957 also
X      missing
ae4   the monthly international airline passenger rates from 1949 to 1960
X      with all observations for the last six years missing
means a matrix that contains the yearly mean values of the logarithmic 
X      series for the years 1955-1960. The data is in the form
X      5.6393   73   84
X      5.7844   85   96
X      ...
X      6.1542  133  144
phi   the phi polynomial of the model. In this case, it is
X      1   1
X      1  12
X     -1  13
X      (see the model above)
theta the theta polynomial of the model. In this case, it is
X      0.396    1
X      0.614   12
X     -0.2431  13
X
The missing observations are given as NaN in the file.
X
The model above was made to logarithmic data.  Therefore, the
following transformed variables should be taken into use:
X
ae1l=log(ae1);
ae2l=log(ae2);
ae3l=log(ae3);
X
Now, the programs can be run. For example, you may call
[zhat,ahat]=inter1(ae1l,phi,theta);
ae2l should be handled with inter2. This is the case also with 
ae3l. Note in the latter case that to get appropriate values
for Julys, you should do the modifications of initial values
described in Karanta (1993), and use the program inter2r.m.
X
To use the program inter1m.m you need input the function with
ae4l, means, phi and theta. The matrix ae4l can be obtained by
the following commands which replace the last six years of the
data with NaN:
X
ae4l=log(aev);
for i=73:144 ae4l(i)=NaN; end;
X
X
Miscellaneous
---------------
X
Documentation on  Matlab may be found in Matlab manuals (MathWorks
1992 a and b).
X
Matlab is a trademark of the MathWorks, Inc.
X
All questions, comments, suggestions for improvement etc. are welcomed
to the Internet e-mail address ilkka.karanta@hut.fi, or
X
Ilkka Karanta
Helsinki Univ. of Technology
Knowledge Engineering Laboratory
Otakaari 1
02150 Espoo
FINLAND
X
X
References
----------
Box and Jenkins: "Time series analysis - forecasting and control",
Holden-Day, San Fransisco 1976.
X
Ilkka Karanta: "Optimization methods in time series interpolation", to
appear in Communications in Statistics, Simulation and Computation.
X
The MathWorks, Inc.: "Matlab Reference Guide", The Mathworks Inc., Natick,
Mass., USA, 1992.
The MathWorks, Inc.: "Matlab User's Guide", The Mathworks Inc., Natick,
Mass., USA, 1992.
X
SHAR_EOF
chmod 0644 README.DOC ||
echo 'restore of README.DOC failed'
Wc_c="`wc -c < 'README.DOC'`"
test 8511 -eq "$Wc_c" ||
	echo 'README.DOC: original size 8511, current size' "$Wc_c"
fi
# ============= aep.m ==============
if test -f 'aep.m' -a X"$1" != X"-c"; then
	echo 'x - skipping aep.m (File already exists)'
else
echo 'x - extracting aep.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'aep.m' &&
ae1=[1.1200000e+02;
X   1.1800000e+02;
X   1.3200000e+02;
X   1.2900000e+02;
X   1.2100000e+02;
X   1.3500000e+02;
X   1.4800000e+02;
X   1.4800000e+02;
X   1.3600000e+02;
X   1.1900000e+02;
X   1.0400000e+02;
X   1.1800000e+02;
X   1.1500000e+02;
X   1.2600000e+02;
X   1.4100000e+02;
X   1.3500000e+02;
X   1.2500000e+02;
X   1.4900000e+02;
X   1.7000000e+02;
X   1.7000000e+02;
X   1.5800000e+02;
X   1.3300000e+02;
X   1.1400000e+02;
X   1.4000000e+02;
X   1.4500000e+02;
X   1.5000000e+02;
X   1.7800000e+02;
X   1.6300000e+02;
X   1.7200000e+02;
X   1.7800000e+02;
X   1.9900000e+02;
X   1.9900000e+02;
X   1.8400000e+02;
X   1.6200000e+02;
X   1.4600000e+02;
X   1.6600000e+02;
X   1.7100000e+02;
X   1.8000000e+02;
X   1.9300000e+02;
X   1.8100000e+02;
X   1.8300000e+02;
X   2.1800000e+02;
X   2.3000000e+02;
X   2.4200000e+02;
X   2.0900000e+02;
X   1.9100000e+02;
X   1.7200000e+02;
X   1.9400000e+02;
X   1.9600000e+02;
X   1.9600000e+02;
X   2.3600000e+02;
X   2.3500000e+02;
X   2.2900000e+02;
X   2.4300000e+02;
X   2.6400000e+02;
X   2.7200000e+02;
X   2.3700000e+02;
X   2.1100000e+02;
X   1.8000000e+02;
X   2.0100000e+02;
X   2.0400000e+02;
X   1.8800000e+02;
X   2.3500000e+02;
X   2.2700000e+02;
X   2.3400000e+02;
X   2.6400000e+02;
X   3.0200000e+02;
X   2.9300000e+02;
X   2.5900000e+02;
X   2.2900000e+02;
X   2.0300000e+02;
X   2.2900000e+02;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X   2.7800000e+02;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X   3.0600000e+02;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X   3.3600000e+02;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X   3.3700000e+02;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X   4.0500000e+02;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X            NaN;
X   4.3200000e+02];
ae2=[1.1200000e+02;
X   1.1800000e+02;
X   1.3200000e+02;
X   1.2900000e+02;
X   1.2100000e+02;
X   1.3500000e+02;
X            NaN;
X   1.4800000e+02;
X   1.3600000e+02;
X   1.1900000e+02;
X   1.0400000e+02;
X   1.1800000e+02;
X   1.1500000e+02;
X   1.2600000e+02;
X   1.4100000e+02;
X   1.3500000e+02;
X   1.2500000e+02;
X   1.4900000e+02;
X   1.7000000e+02;
X   1.7000000e+02;
X   1.5800000e+02;
X   1.3300000e+02;
X   1.1400000e+02;
X   1.4000000e+02;
X   1.4500000e+02;
X   1.5000000e+02;
X   1.7800000e+02;
X   1.6300000e+02;
X   1.7200000e+02;
X   1.7800000e+02;
X   1.9900000e+02;
X   1.9900000e+02;
X   1.8400000e+02;
X   1.6200000e+02;
X   1.4600000e+02;
X   1.6600000e+02;
X   1.7100000e+02;
X   1.8000000e+02;
X   1.9300000e+02;
X   1.8100000e+02;
X   1.8300000e+02;
X   2.1800000e+02;
X   2.3000000e+02;
X   2.4200000e+02;
X   2.0900000e+02;
X   1.9100000e+02;
X   1.7200000e+02;
X   1.9400000e+02;
X   1.9600000e+02;
X   1.9600000e+02;
X   2.3600000e+02;
X   2.3500000e+02;
X   2.2900000e+02;
X   2.4300000e+02;
X   2.6400000e+02;
X   2.7200000e+02;
X   2.3700000e+02;
X   2.1100000e+02;
X   1.8000000e+02;
X   2.0100000e+02;
X   2.0400000e+02;
X   1.8800000e+02;
X   2.3500000e+02;
X   2.2700000e+02;
X   2.3400000e+02;
X   2.6400000e+02;
X   3.0200000e+02;
X   2.9300000e+02;
X   2.5900000e+02;
X   2.2900000e+02;
X   2.0300000e+02;
X   2.2900000e+02;
X   2.4200000e+02;
X   2.3300000e+02;
X   2.6700000e+02;
X   2.6900000e+02;
X   2.7000000e+02;
X   3.1500000e+02;
X   3.6400000e+02;
X   3.4700000e+02;
X   3.1200000e+02;
X   2.7400000e+02;
X   2.3700000e+02;
X   2.7800000e+02;
X   2.8400000e+02;
X   2.7700000e+02;
X   3.1700000e+02;
X   3.1300000e+02;
X   3.1800000e+02;
X   3.7400000e+02;
X   4.1300000e+02;
X   4.0500000e+02;
X   3.5500000e+02;
X   3.0600000e+02;
X   2.7100000e+02;
X   3.0600000e+02;
X   3.1500000e+02;
X   3.0100000e+02;
X   3.5600000e+02;
X   3.4800000e+02;
X   3.5500000e+02;
X            NaN;
X            NaN;
X            NaN;
X   4.0400000e+02;
X   3.4700000e+02;
X   3.0500000e+02;
X   3.3600000e+02;
X   3.4000000e+02;
X   3.1800000e+02;
X   3.6200000e+02;
X   3.4800000e+02;
X   3.6300000e+02;
X   4.3500000e+02;
X   4.9100000e+02;
X   5.0500000e+02;
X   4.0400000e+02;
X   3.5900000e+02;
X   3.1000000e+02;
X   3.3700000e+02;
X   3.6000000e+02;
X   3.4200000e+02;
X   4.0600000e+02;
X   3.9600000e+02;
X   4.2000000e+02;
X   4.7200000e+02;
X   5.4800000e+02;
X   5.5900000e+02;
X   4.6300000e+02;
X   4.0700000e+02;
X   3.6200000e+02;
X   4.0500000e+02;
X   4.1700000e+02;
X   3.9100000e+02;
X   4.1900000e+02;
X   4.6100000e+02;
X   4.7200000e+02;
X   5.3500000e+02;
X            NaN;
X   6.0600000e+02;
X   5.0800000e+02;
X   4.6100000e+02;
X   3.9000000e+02;
X   4.3200000e+02];
ae3=[1.1200000e+02;
X   1.1800000e+02;
X   1.3200000e+02;
X   1.2900000e+02;
X   1.2100000e+02;
X   1.3500000e+02;
X            NaN;
X   1.4800000e+02;
X   1.3600000e+02;
X   1.1900000e+02;
X   1.0400000e+02;
X   1.1800000e+02;
X   1.1500000e+02;
X   1.2600000e+02;
X   1.4100000e+02;
X   1.3500000e+02;
X   1.2500000e+02;
X   1.4900000e+02;
X            NaN;
X   1.7000000e+02;
X   1.5800000e+02;
X   1.3300000e+02;
X   1.1400000e+02;
X   1.4000000e+02;
X   1.4500000e+02;
X   1.5000000e+02;
X   1.7800000e+02;
X   1.6300000e+02;
X   1.7200000e+02;
X   1.7800000e+02;
X            NaN;
X   1.9900000e+02;
X   1.8400000e+02;
X   1.6200000e+02;
X   1.4600000e+02;
X   1.6600000e+02;
X   1.7100000e+02;
X   1.8000000e+02;
X   1.9300000e+02;
X   1.8100000e+02;
X   1.8300000e+02;
X   2.1800000e+02;
X            NaN;
X   2.4200000e+02;
X   2.0900000e+02;
X   1.9100000e+02;
X   1.7200000e+02;
X   1.9400000e+02;
X   1.9600000e+02;
X   1.9600000e+02;
X   2.3600000e+02;
X   2.3500000e+02;
X   2.2900000e+02;
X   2.4300000e+02;
X            NaN;
X   2.7200000e+02;
X   2.3700000e+02;
X   2.1100000e+02;
X   1.8000000e+02;
X   2.0100000e+02;
X   2.0400000e+02;
X   1.8800000e+02;
X   2.3500000e+02;
X   2.2700000e+02;
X   2.3400000e+02;
X   2.6400000e+02;
X            NaN;
X   2.9300000e+02;
X   2.5900000e+02;
X   2.2900000e+02;
X   2.0300000e+02;
X   2.2900000e+02;
X   2.4200000e+02;
X   2.3300000e+02;
X   2.6700000e+02;
X   2.6900000e+02;
X   2.7000000e+02;
X   3.1500000e+02;
X            NaN;
X   3.4700000e+02;
X   3.1200000e+02;
X   2.7400000e+02;
X   2.3700000e+02;
X   2.7800000e+02;
X   2.8400000e+02;
X   2.7700000e+02;
X   3.1700000e+02;
X   3.1300000e+02;
X   3.1800000e+02;
X   3.7400000e+02;
X            NaN;
X   4.0500000e+02;
X   3.5500000e+02;
X   3.0600000e+02;
X   2.7100000e+02;
X   3.0600000e+02;
X   3.1500000e+02;
X   3.0100000e+02;
X   3.5600000e+02;
X   3.4800000e+02;
X   3.5500000e+02;
X            NaN;
X            NaN;
X            NaN;
X   4.0400000e+02;
X   3.4700000e+02;
X   3.0500000e+02;
X   3.3600000e+02;
X   3.4000000e+02;
X   3.1800000e+02;
X   3.6200000e+02;
X   3.4800000e+02;
X   3.6300000e+02;
X   4.3500000e+02;
X            NaN;
X   5.0500000e+02;
X   4.0400000e+02;
X   3.5900000e+02;
X   3.1000000e+02;
X   3.3700000e+02;
X   3.6000000e+02;
X   3.4200000e+02;
X   4.0600000e+02;
X   3.9600000e+02;
X   4.2000000e+02;
X   4.7200000e+02;
X            NaN;
X   5.5900000e+02;
X   4.6300000e+02;
X   4.0700000e+02;
X   3.6200000e+02;
X   4.0500000e+02;
X   4.1700000e+02;
X   3.9100000e+02;
X   4.1900000e+02;
X   4.6100000e+02;
X   4.7200000e+02;
X   5.3500000e+02;
X            NaN;
X   6.0600000e+02;
X   5.0800000e+02;
X   4.6100000e+02;
X   3.9000000e+02;
X   4.3200000e+02];
aev=[1.1200000e+02;
X   1.1800000e+02;
X   1.3200000e+02;
X   1.2900000e+02;
X   1.2100000e+02;
X   1.3500000e+02;
X   1.4800000e+02;
X   1.4800000e+02;
X   1.3600000e+02;
X   1.1900000e+02;
X   1.0400000e+02;
X   1.1800000e+02;
X   1.1500000e+02;
X   1.2600000e+02;
X   1.4100000e+02;
X   1.3500000e+02;
X   1.2500000e+02;
X   1.4900000e+02;
X   1.7000000e+02;
X   1.7000000e+02;
X   1.5800000e+02;
X   1.3300000e+02;
X   1.1400000e+02;
X   1.4000000e+02;
X   1.4500000e+02;
X   1.5000000e+02;
X   1.7800000e+02;
X   1.6300000e+02;
X   1.7200000e+02;
X   1.7800000e+02;
X   1.9900000e+02;
X   1.9900000e+02;
X   1.8400000e+02;
X   1.6200000e+02;
X   1.4600000e+02;
X   1.6600000e+02;
X   1.7100000e+02;
X   1.8000000e+02;
X   1.9300000e+02;
X   1.8100000e+02;
X   1.8300000e+02;
X   2.1800000e+02;
X   2.3000000e+02;
X   2.4200000e+02;
X   2.0900000e+02;
X   1.9100000e+02;
X   1.7200000e+02;
X   1.9400000e+02;
X   1.9600000e+02;
X   1.9600000e+02;
X   2.3600000e+02;
X   2.3500000e+02;
X   2.2900000e+02;
X   2.4300000e+02;
X   2.6400000e+02;
X   2.7200000e+02;
X   2.3700000e+02;
X   2.1100000e+02;
X   1.8000000e+02;
X   2.0100000e+02;
X   2.0400000e+02;
X   1.8800000e+02;
X   2.3500000e+02;
X   2.2700000e+02;
X   2.3400000e+02;
X   2.6400000e+02;
X   3.0200000e+02;
X   2.9300000e+02;
X   2.5900000e+02;
X   2.2900000e+02;
X   2.0300000e+02;
X   2.2900000e+02;
X   2.4200000e+02;
X   2.3300000e+02;
X   2.6700000e+02;
X   2.6900000e+02;
X   2.7000000e+02;
X   3.1500000e+02;
X   3.6400000e+02;
X   3.4700000e+02;
X   3.1200000e+02;
X   2.7400000e+02;
X   2.3700000e+02;
X   2.7800000e+02;
X   2.8400000e+02;
X   2.7700000e+02;
X   3.1700000e+02;
X   3.1300000e+02;
X   3.1800000e+02;
X   3.7400000e+02;
X   4.1300000e+02;
X   4.0500000e+02;
X   3.5500000e+02;
X   3.0600000e+02;
X   2.7100000e+02;
X   3.0600000e+02;
X   3.1500000e+02;
X   3.0100000e+02;
X   3.5600000e+02;
X   3.4800000e+02;
X   3.5500000e+02;
X   4.2200000e+02;
X   4.6500000e+02;
X   4.6700000e+02;
X   4.0400000e+02;
X   3.4700000e+02;
X   3.0500000e+02;
X   3.3600000e+02;
X   3.4000000e+02;
X   3.1800000e+02;
X   3.6200000e+02;
X   3.4800000e+02;
X   3.6300000e+02;
X   4.3500000e+02;
X   4.9100000e+02;
X   5.0500000e+02;
X   4.0400000e+02;
X   3.5900000e+02;
X   3.1000000e+02;
X   3.3700000e+02;
X   3.6000000e+02;
X   3.4200000e+02;
X   4.0600000e+02;
X   3.9600000e+02;
X   4.2000000e+02;
X   4.7200000e+02;
X   5.4800000e+02;
X   5.5900000e+02;
X   4.6300000e+02;
X   4.0700000e+02;
X   3.6200000e+02;
X   4.0500000e+02;
X   4.1700000e+02;
X   3.9100000e+02;
X   4.1900000e+02;
X   4.6100000e+02;
X   4.7200000e+02;
X   5.3500000e+02;
X   6.2200000e+02;
X   6.0600000e+02;
X   5.0800000e+02;
X   4.6100000e+02;
X   3.9000000e+02;
X   4.3200000e+02];
means=[5.6392685e+00   7.3000000e+01   8.4000000e+01;
X   5.7844288e+00   8.5000000e+01   9.6000000e+01;
X   5.8983724e+00   9.7000000e+01   1.0800000e+02;
X   5.9305823e+00   1.0900000e+02   1.2000000e+02;
X   6.0483677e+00   1.2100000e+02   1.3200000e+02;
X   6.1542150e+00   1.3300000e+02   1.4400000e+02];
phi=[1.0000000e+00   1.0000000e+00;
X   1.0000000e+00   1.2000000e+01;
X  -1.0000000e+00   1.3000000e+01];
theta=[3.9600000e-01   1.0000000e+00;
X   6.1400000e-01   1.2000000e+01;
X  -2.4314400e-01   1.3000000e+01];
SHAR_EOF
chmod 0666 aep.m ||
echo 'restore of aep.m failed'
Wc_c="`wc -c < 'aep.m'`"
test 10812 -eq "$Wc_c" ||
	echo 'aep.m: original size 10812, current size' "$Wc_c"
fi
# ============= app.m ==============
if test -f 'app.m' -a X"$1" != X"-c"; then
	echo 'x - skipping app.m (File already exists)'
else
echo 'x - extracting app.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'app.m' &&
function x = app(phi,direct,xvec,place);
X
% Calculates the value of phi(B)z_t in the lag polynomial
% (1-phi(B))z_t=(1-phi_1B-phi2B^2-...)z_t
% The call is of the form
%   x = app(phi,direct,xvec,place);
% where phi includes the parameter values of the polynomial phi 
%         at the places phi(:,1)
%       the lags corresponding to the parameter values at the places phi(:,2)
%       direct is the direction of application: -1 = B, 1 = F.
%         Here Bz_t=z_{t-1}, and Fz_t=z_{t+1}.
%       xvec is the time series to be filtered
%       place is the value of the time variable, that is, it is 
%         the place in the time series of current interest
X
% Ilkka Karanta 1991-1993
X
phisiz=size(phi);
parnum=phisiz(1);
xlen=length(xvec);
X
x=0.0;
if (length(phi)>0)
for i=1:parnum
X if (place+direct*phi(i,2)<=xlen) & (place+direct*phi(i,2)>=1)
X   x=x+phi(i,1)*xvec(place+direct*phi(i,2));
X end
end
end
X
SHAR_EOF
chmod 0644 app.m ||
echo 'restore of app.m failed'
Wc_c="`wc -c < 'app.m'`"
test 908 -eq "$Wc_c" ||
	echo 'app.m: original size 908, current size' "$Wc_c"
fi
# ============= autocorr.m ==============
if test -f 'autocorr.m' -a X"$1" != X"-c"; then
	echo 'x - skipping autocorr.m (File already exists)'
else
echo 'x - extracting autocorr.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'autocorr.m' &&
function [p,psi]=autocorr(phi,theta,n);
X
% AUTOCORR calculates the autocorrelations of the ARIMA model
% phi(B)z_t=theta(B)a_t. The call is of the form
%  [p psi]=autocorr(phi,theta,n);
% where p     is a vector containing the autocorrelations
%       psi   is theta(B)/phi(B)
%       phi   is the AR polynomial of the model
%       theta is the MA polynomial of the model
%       n     is the number of lags wished
% NOTE  length(pi)=n+1, since the first element of p is the 
% zeroth lag autocorrelation
X
% Reference: P.J. Brockwell and R.A. Davis: "Time series: 
% theory and methods"  (Springer 1987), sec. 3.3
X
% Ilkka Karanta 1991-1993
X
apuphi=zeros(n,1);
aputheta=zeros(n,1);
phisize=size(phi);
thetasize=size(theta);
for i=1:phisize(1)
X  apuphi(phi(i,2))=phi(i,1);
end
for i=1:thetasize(1)
X  aputheta(theta(i,2))=theta(i,1);
end
X
if phisize(1)>0
X  philen=phi(phisize(1),2);
else
X  philen=0;
end
if thetasize(1)>0
X  thetalen=theta(thetasize(1),2);
else
X  thetalen = 0;
end
maxlen=max(philen,thetalen+1)+1
X
% calculation of the psi vector
psi=pttopi(theta,phi,n);
X
% calculation of the p vector
% first the initial values
kerrmat=eye(maxlen);
for i=1:maxlen
X  for j=1:i-1
X    kerrmat(i,j)=kerrmat(i,j)-apuphi(i-j);
X  end
X  for j=2:maxlen
X    kerrmat(i,j)=kerrmat(i,j)-apuphi(i+j-2);
X  end
end
kerrmat
rhsvek=zeros(maxlen,1);
rhsvek(1)=1;
for j=1:thetalen
X  rhsvek(1)=rhsvek(1)-aputheta(j)*psi(j);
end
for i=1:maxlen-1
X  rhsvek(1+i)=-aputheta(i);
X  for j=i+1:thetalen
X    rhsvek(1+i)=rhsvek(1+i)-aputheta(j)*psi(j-i);
X  end
end
rhsvek
X
p=zeros(n+1,1);
p(1:maxlen)=kerrmat\rhsvek;
X
% the rest of the difference equation
X
for t=maxlen+1:n+1
X  p(t)=app(phi,-1,p,t);
end
SHAR_EOF
chmod 0644 autocorr.m ||
echo 'restore of autocorr.m failed'
Wc_c="`wc -c < 'autocorr.m'`"
test 1672 -eq "$Wc_c" ||
	echo 'autocorr.m: original size 1672, current size' "$Wc_c"
fi
# ============= broyden.m ==============
if test -f 'broyden.m' -a X"$1" != X"-c"; then
	echo 'x - skipping broyden.m (File already exists)'
else
echo 'x - extracting broyden.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'broyden.m' &&
function H = broyden(pnew,p,H,graduu,gradvan);
% calculate the new approximation of the Hessian matrix by Broyden's method.
% the call is of the form
%    H = broyden(pnew,p,H,graduu,gradvan);
% Here pnew    is the vector of new parameter values, 
%      p       is the vector of old parameter values,
%      H       is the Hessian matrix of the cost function,
%      graduu  is the gradient of the cost function at pnew,
%      gradvan is the gradient of the cost function at p
X
% Ilkka Karanta 1991-1993
X
pdif=pnew-p;
q=graduu-gradvan;
apu=pdif-H*q;
% rank one correction
H=H+apu*apu'/(q'*apu);
SHAR_EOF
chmod 0644 broyden.m ||
echo 'restore of broyden.m failed'
Wc_c="`wc -c < 'broyden.m'`"
test 597 -eq "$Wc_c" ||
	echo 'broyden.m: original size 597, current size' "$Wc_c"
fi
# ============= forec.m ==============
if test -f 'forec.m' -a X"$1" != X"-c"; then
	echo 'x - skipping forec.m (File already exists)'
else
echo 'x - extracting forec.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'forec.m' &&
function zhat=forec(z,a,phi,theta,flength)
X
% FOREC calculates forecasts for time series z up to lead time
% flength. An ARIMA model is assumed to be known to the series.
% The call is of the form
%    zhat=forec(z,a,phi,theta,flength)
% where
%   z     the history of the time series to be forecasted
%   a     estimates of noise corresponding to z
%   phi   the AR polynomial of the model
%   theta the MA polynomial of the model
%   flength forecasting lead time
%   zhat    time series containing both z and the forecasts
X
% Ilkka Karanta 1991-1993
X
tyhjat=max(phi(length(phi),2), theta(length(theta),2));
slen=length(z);
zhat=zeros(slen+flength,1);
zhat(1:slen)=z;
ahat=zeros(slen+flength,1);
X
X
% initialization of noises
if a'*a==0  
for j=tyhjat+1:slen
%   zhat(j)=app(phi,-1,z,j)-app(theta,-1,ahat,j);
X   ahat(j)=z(j)-app(phi,-1,z,j)+app(theta,-1,ahat,j);
end
for i=slen-tyhjat-1:-1:1
%    zhat(i)=app(phi,1,z,i)-app(theta,1,ahat,i);
X    ahat(i)=z(i)-app(phi,1,z,i)+app(theta,1,ahat,i);
end
for j=tyhjat+1:slen
%      zhat(j)=app(phi,-1,z,j)-app(theta,-1,ahat,j);
X      ahat(j)=z(j)-app(phi,-1,z,j)+app(theta,-1,ahat,j);
end
else
X  ahat(1:slen)=a;
end
X
% the actual forecasting
X
for t=slen+1:slen+flength
X   zhat(t)=app(phi,-1,zhat,t)-app(theta,-1,ahat,t);
end
SHAR_EOF
chmod 0644 forec.m ||
echo 'restore of forec.m failed'
Wc_c="`wc -c < 'forec.m'`"
test 1268 -eq "$Wc_c" ||
	echo 'forec.m: original size 1268, current size' "$Wc_c"
fi
# ============= fval.m ==============
if test -f 'fval.m' -a X"$1" != X"-c"; then
	echo 'x - skipping fval.m (File already exists)'
else
echo 'x - extracting fval.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fval.m' &&
function [zhat,ahat,f]=fval(aest,z,zhat,m,n,tyhjat,phi,theta)
% FVAL calculates the estimated z, a and f (the norm of a)
% with known estimates of noise aest. The call is of the form
%  [zhat,ahat,f]=fval(aest,z,zhat,m,n,tyhjat,phi,theta)
% where
%   zhat includes estimates of observations in the time series (the
%        same as the actual observation when known)
%   ahat the corresponding noises
%   aest noise estimates for the missing values of the time series
%   z    the time series the missing values of which are being calculated
%   m    length of the portion of the time series at the beginning
%        with no missing observations
%   n    length of the rest of the time series
%   tyhjat max(order of phi, order of theta)
%   phi  the phi parameters of the time series
%   theta the theta parameters of the time series
X
% Ilkka Karanta 1991-1993
X
ahat=zeros(size(zhat));
apun=find(isnan(z));
ahat(1:tyhjat)=zeros(tyhjat,1);
X for i=1:length(apun)
X    ahat(apun(i))=aest(i);
X end
X for j=tyhjat+1:n+m
X   if isnan(z(j))
X      zhat(j)=app(phi,-1,zhat,j)+ahat(j)-app(theta,-1,ahat,j);
X   else
X      ahat(j)=zhat(j)-app(phi,-1,zhat,j)+app(theta,-1,ahat,j);
X   end
X end
X  for i=m+n-tyhjat-1:-1:1
X    if isnan(z(i))
X      zhat(i)=app(phi,1,zhat,i)+ahat(i)-app(theta,1,ahat,i);
X    else
X      ahat(i)=zhat(i)-app(phi,1,zhat,i)+app(theta,1,ahat,i);
X    end;
X  end
X for j=tyhjat+1:n+m
X   if isnan(z(j))
X      zhat(j)=app(phi,-1,zhat,j)+ahat(j)-app(theta,-1,ahat,j);
X   else
X      ahat(j)=zhat(j)-app(phi,-1,zhat,j)+app(theta,-1,ahat,j);
X   end
X end
X
X f=norm(ahat);
SHAR_EOF
chmod 0644 fval.m ||
echo 'restore of fval.m failed'
Wc_c="`wc -c < 'fval.m'`"
test 1569 -eq "$Wc_c" ||
	echo 'fval.m: original size 1569, current size' "$Wc_c"
fi
# ============= inter1.m ==============
if test -f 'inter1.m' -a X"$1" != X"-c"; then
	echo 'x - skipping inter1.m (File already exists)'
else
echo 'x - extracting inter1.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'inter1.m' &&
function [zhat,ahat]=inter1(z,phi,theta)
X
% INTER1 calculates estimates for missing observations
% when an ARIMA model 
%    phi(B)z_t=theta(B)a_t
% for the time series is known.
% This version of the algorithm is specially adapted to handling time
% series of the form
%   [Z U1 z1 U2 z2 U3 ... zr],
% where Z is a sequence of known observations (it might contain
% no elements); U1, U2, ..., Ur are sequences of unknown observations;
% and z1, z2, ..., zr are sequences of one (or at most a few) observations.
% The call is of the form
%   [zhat,ahat]=inter1(z,phi,theta)
% where
%   z     the time series that contains missing observations.
%         The missing observations must be supplied as NaN.
%         An example series: [1 3.1 5 4 NaN 12.7 NaN NaN 3 NaN 1]
%   phi   the AR polynomial of the model; It should be in the form
%           [phi_1 lag_1; phi_2 lag_2; ... phi_p lag_p],
%         where phi_i's are the parameter values and lag_i's are the 
%         corresponding lags in the model.
%   theta the MA polynomial of the model. The format should be as with phi.
%   zhat   time series containing both the known observations and
%          estimates for the unknown ones
%   ahat   the corresponding noise values
X
% The function is an implementation of Algorithm 1 in the paper
%   Ilkka Karanta: "Optimization methods in time series interpolation",
%   Communications in Statistics, Simulation and Computation (to appear).
X
% Ilkka Karanta 1991-1993
X
% internal variables
%   lambda the Lagrangean multipliers corresponding to the ARIMA equation
%   p      the Lagrangean multipliers corresponding to known observations
%   grad   the gradient of the cost function with respect to p
%   n      length of the part of time series from the first missing
%          observations to end
%   m      length of the part of time series from the beginning to
%          the observation preceeding the first unknown observation
%   count  the number of iteration counts performed so far. Can be
%          used for performance screening purposes.
%
X
ukind=find(isnan(z));
m=ukind(1)-1;
n=length(z)-m;
initzer=max(phi(length(phi),2),theta(length(theta),2));
zhat=zeros(n+m,1);
zhat(1:m)=z(1:m);
ahat=zeros(n+m,1);
lambda=zeros(n+m,1);
apu=z(find(~isnan(z(m+1:m+n))));
H=eye(length(apu));
p=zeros(length(apu),1);
gradnew=zeros(length(apu),1);
gradold=zeros(length(apu),1);
X
% initialization of the noise variables
for t=initzer+1:m
X   ahat(t)=zhat(t)-app(phi,-1,zhat,t)+app(theta,-1,ahat,t);
end
for t=m-initzer-1:1
X   ahat(t)=zhat(t)-app(phi,1,zhat,t)+app(theta,1,ahat,t);
end
for t=initzer+1:m
X   ahat(t)=zhat(t)-app(phi,-1,zhat,t)+app(theta,-1,ahat,t);
end
X
k=0;
for j=1:n
X   zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   if ~isnan(z(m+j))
X      k=k+1;
X      gradnew(k)=z(m+j)-zhat(m+j);
X   end
end
gradold=0.9*gradnew;
epsilon=gradnew'*gradnew;
X
count=0;
while epsilon>1.0e-12
X count=count+1;
X pvan=p;
X p=p-H*gradnew;
X k=length(apu);
X for i=n:-1:1
X   if isnan(z(m+i)) 
X      lambda(m+i)=app(phi,1,lambda,m+i);
X   else
X      lambda(m+i)=app(phi,1,lambda,m+i)+p(k);
X      k=k-1;
X   end
X   ahat(m+i)=lambda(m+i)-app(theta,1,lambda,m+i);
X end
X k=0;
X gradold=gradnew;
X for j=1:n
X   zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   if ~isnan(z(m+j))
X      k=k+1;
X      gradnew(k)=z(m+j)-zhat(m+j);
X   end
X end
X H=broyden(p,pvan,H,gradnew,gradold);
X epsilon=gradnew'*gradnew;
end
X
X
X
SHAR_EOF
chmod 0644 inter1.m ||
echo 'restore of inter1.m failed'
Wc_c="`wc -c < 'inter1.m'`"
test 3429 -eq "$Wc_c" ||
	echo 'inter1.m: original size 3429, current size' "$Wc_c"
fi
# ============= inter1m.m ==============
if test -f 'inter1m.m' -a X"$1" != X"-c"; then
	echo 'x - skipping inter1m.m (File already exists)'
else
echo 'x - extracting inter1m.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'inter1m.m' &&
function [zhat,ahat]=inter1m(z,a,means,phi,theta)
X
% count estimates of missing values in a time series. Means of 
% the missing values for each interval of missing values is 
% assumed to be known.
%
% input variables:
%   z     time series with missing values
%   a     noise values of z with respect to the model phi(B)z_t=theta(B)a_t 
%   means a matrix that includes 
%     1. the mean of the time series on the given interval
%     2. the starting point of the given interval
%     3. the end point of the given interval
%   phi   the AR polynomial of the model
%   theta the MA polynomial of the model
%
% output variables:
%   zhat  a time series that includes the known values of z and estimates
%         for missing values
%   ahat  noise values corresponding to zhat
X
% Method see the paper Ilkka Karanta: "Optimization methods in
% time series missing values interpolation". 
X
% Ilkka Karanta 1991-1993
X
% internal variables
%   lambda co-state variable (Lagrangean multiplier) corresponding 
%          to the model equation
%   p      co-state variable corresponding to each interval of missing 
%          values
%   grad   gradient of the cost function with respect to p
%   n      length of the part of time series with missing values
%   m      length of the part of time series with no missing values
X
ukind=find(isnan(z));
m=ukind(1)-1;
n=length(z)-m;
mkoko=size(means);
tyhjat=max(theta(length(theta),2),phi(length(phi),2));
zhat=zeros(n+m,1);
zhat(1:m)=z(1:m);
ahat=zeros(n+m,1);
ahat(1:m)=a(1:m);
lambda=zeros(n+m,1);
H=eye(mkoko(1));
p=zeros(mkoko(1),1);
graduu=zeros(mkoko(1),1);
gradvan=zeros(mkoko(1),1);
X
k=0;
for j=1:n
X   zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   if ~isnan(z(m+j))
X      ahat(m+j)=z(m+j)-zhat(m+j);
X      zhat(m+j)=z(m+j);
X   end
end
for i=1:mkoko(1)
X   graduu(i)=means(i,1)-mean(zhat(means(i,2):means(i,3)));
end
gradvan=0.9*graduu;
graduu
epsilon=graduu'*graduu
p
H
X
count=0;
while epsilon>1.0e-12
X count=count+1;
X pvan=p;
X p=p-H*graduu;
X t=m+n;
X i=mkoko(1);
X while t>1
X   if t==means(i,3)
X      for t=means(i,3):-1:means(i,2)
X         lambda(t)=app(phi,1,lambda,t)+p(i)/(means(i,3)-means(i,2)+1);
X         ahat(t)=lambda(t)-app(theta,1,lambda,t);
X      end
% t==means(i,2)
X      t=t-1;
X      if i>1
X         i=i-1;
X      end
X   else
X      lambda(t)=app(phi,1,lambda,t);
X      ahat(t)=lambda(t)-app(theta,1,lambda,t);
X      t=t-1;
X   end
X end
X gradvan=graduu;
X for t=tyhjat+1:m+n
X   if isnan(z(t))
X      zhat(t)=app(phi,-1,zhat,t)+ahat(t)-app(theta,-1,ahat,t);
X   else
X      ahat(t)=z(t)-app(phi,-1,zhat,t)+app(theta,-1,ahat,t);
X   end
X end
X for i=1:mkoko(1)
X    graduu(i)=means(i,1)-mean(zhat(means(i,2):means(i,3)));
X end
X
X H=broyden(p,pvan,H,graduu,gradvan);
X epsilon=graduu'*graduu;
end
X
SHAR_EOF
chmod 0644 inter1m.m ||
echo 'restore of inter1m.m failed'
Wc_c="`wc -c < 'inter1m.m'`"
test 2767 -eq "$Wc_c" ||
	echo 'inter1m.m: original size 2767, current size' "$Wc_c"
fi
# ============= inter1r.m ==============
if test -f 'inter1r.m' -a X"$1" != X"-c"; then
	echo 'x - skipping inter1r.m (File already exists)'
else
echo 'x - extracting inter1r.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'inter1r.m' &&
function [zhat,ahat]=inter1r(z,a,zest,phi,theta,tol)
X
% INTER1R calculates estimates for missing observations
% when an ARIMA model 
%    phi(B)z_t=theta(B)a_t
% for the time series is known. It is, actually, inter1 with more 
% input variables. This version of the algorithm is specially adapted 
% to handling time series of the form
%   [Z U1 z1 U2 z2 U3 ... zr],
% where Z is a sequence of known observations (it might contain
% no elements); U1, U2, ..., Ur are sequences of unknown observations;
% and z1, z2, ..., zr are sequences of one (or at most a few) observations.
% The call is of the form
%   [zhat,ahat]=inter1r(z,a,zest,phi,theta,tol)
% where
%   z     the time series that contains missing observations.
%         The missing observations must be supplied as NaN.
%         An example series: [1 3.1 5 4 NaN 12.7 NaN NaN 3 NaN 1]
%   a     initial estimates for the noise variables. May be set to zero.
%   zest  the known time series values and initial estimates for the missing
%         observations
%   phi   the AR polynomial of the model; It should be in the form
%           [phi_1 lag_1; phi_2 lag_2; ... phi_p lag_p],
%         where phi_i's are the parameter values and lag_i's are the 
%         corresponding lags in the model.
%   theta the MA polynomial of the model. The format should be as with phi.
%   tol   tolerance in the stopping criterion; the maximum
%         allowed square of the norm of gradient
%   zhat   time series containing both the known observations and
%          estimates for the unknown ones
%   ahat   the corresponding noise values
X
% The function is an implementation of Algorithm 1 in the paper
%   Ilkka Karanta: "Optimization methods in time series interpolation",
%   Communications in Statistics, Simulation and Computation (to appear).
X
% Ilkka Karanta 1993
X
% internal variables
%   lambda the Lagrangean multipliers corresponding to the ARIMA equation
%   p      the Lagrangean multipliers corresponding to known observations
%   grad   the gradient of the cost function with respect to p
%   n      length of the part of time series from the first missing
%          observations to end
%   m      length of the part of time series from the beginning to
%          the observation preceeding the first unknown observation
%
X
ukind=find(isnan(z));
m=ukind(1)-1;
n=length(z)-m;
initzer=max(phi(length(phi),2),theta(length(theta),2));
zhat=zeros(n+m,1);
zhat=z;
zhat(ukind)=zest(ukind);
ahat=zeros(n+m,1);
ahat=a;
lambda=zeros(n+m,1);
apu=z(find(~isnan(z(m+1:m+n))));
H=eye(length(apu));
p=zeros(length(apu),1);
gradnew=zeros(length(apu),1);
gradold=zeros(length(apu),1);
X
% initialization of the noise variables
for t=initzer+1:m
X   ahat(t)=zhat(t)-app(phi,-1,zhat,t)+app(theta,-1,ahat,t);
end
for t=m-initzer-1:1
X   ahat(t)=zhat(t)-app(phi,1,zhat,t)+app(theta,1,ahat,t);
end
for t=initzer+1:m
X   ahat(t)=zhat(t)-app(phi,-1,zhat,t)+app(theta,-1,ahat,t);
end
X
k=0;
for j=1:n
X   zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   if ~isnan(z(m+j))
X      k=k+1;
X      gradnew(k)=z(m+j)-zhat(m+j);
X   end
end
gradold=0.9*gradnew;
epsilon=gradnew'*gradnew
p
H
X
count=0;
while epsilon>tol
X count=count+1;
X pvan=p;
X p=p-H*gradnew;
X k=length(apu);
X for i=n:-1:1
X   if isnan(z(m+i)) 
X      lambda(m+i)=app(phi,1,lambda,m+i);
X   else
X      lambda(m+i)=app(phi,1,lambda,m+i)+p(k);
X      k=k-1;
X   end
X   ahat(m+i)=lambda(m+i)-app(theta,1,lambda,m+i);
X end
X k=0;
X gradold=gradnew;
X for j=1:n
X   zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   if ~isnan(z(m+j))
X      k=k+1;
X      gradnew(k)=z(m+j)-zhat(m+j);
X   end
X end
X H=broyden(p,pvan,H,gradnew,gradold);
X epsilon=gradnew'*gradnew;
end
X
count
X
X
SHAR_EOF
chmod 0644 inter1r.m ||
echo 'restore of inter1r.m failed'
Wc_c="`wc -c < 'inter1r.m'`"
test 3692 -eq "$Wc_c" ||
	echo 'inter1r.m: original size 3692, current size' "$Wc_c"
fi
# ============= inter2.m ==============
if test -f 'inter2.m' -a X"$1" != X"-c"; then
	echo 'x - skipping inter2.m (File already exists)'
else
echo 'x - extracting inter2.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'inter2.m' &&
function [zhat,ahat]=inter2(z,phi,theta)
X
% calculate estimates for the missing values in a time series.
% It is assumed that there are relatively few unknown points.
% The call is of the form
%
%    [zhat,ahat]=inter2(z,phi,theta)
%
% where
%   z     a time series with missing values (see the accompanying README file)
%   phi   the AR polynomial of the model
%   theta the MA polynomial of the model
%   zhat  a time series containing both the known time series values
%         and estimates for the unknown values
%   ahat  the noise estimates corresponding to zhat
X
% For the method see algorithm 2 in the paper 
%   Ilkka Karanta: Optimization methods in time series interpolation
% (Communications in Statistics, to appear)
X
% Ilkka Karanta 1993
X
% internal variables
%   lambda the co-state variable corresponding to the ARIMA equation
%   grad   the gradient of the cost function with respect to the
%          unknown noises a_t
%   n      length of the part of the time series with missing values
%   m      length of the part of the time series in the beginning of
%          the series with no missing values
%
X
ukind=find(isnan(z));
m=ukind(1)-1;
n=length(z)-m;
zhat=zeros(n+m,1);
zhat(1:m)=z(1:m);
ahat=zeros(n+m,1);
% the maximunm length of lambda is 
%     n+m+max(the order of phi, the order of theta)
lambda=zeros(n+m+336,1);
apu=z(find(isnan(z(m+1:m+n))));
H=eye(length(apu));
aest=zeros(length(apu),1);
graduu=zeros(length(apu),1);
gradvan=zeros(length(apu),1);
X
% initialization
k=0;
for j=1:n
X   zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   if isnan(z(m+j))
X      k=k+1;
X      aest(k)=0;
X   else
X      ahat(m+j)=z(m+j)-zhat(m+j);
X      zhat(m+j)=z(m+j);
X   end
end
% k=the number of points with unknown data
for i=n:-1:1
X   if isnan(z(m+i))
X      lambda(m+i)=app(phi,1,lambda,m+i);
X      graduu(k)=ahat(m+i)+app(theta,1,lambda,m+i)-lambda(m+i);
X      k=k-1;
X   else
X      lambda(m+i)=ahat(m+i)+app(theta,1,lambda,m+i);
X   end;
end;
sizeuu=size(graduu);
epsilon=graduu'*graduu;
X
% the actual iteration loop
count=0;
len=length(apu);
while epsilon>1.0e-8
X count=count+1;
X aestvan=aest;
X aest=aest-H*graduu;
X k=0;
X for j=1:n
X   if isnan(z(m+j))
X      k=k+1;
X      ahat(m+j)=aest(k);
X      zhat(m+j)=app(phi,-1,zhat,m+j)+ahat(m+j)-app(theta,-1,ahat,m+j);
X   else
X      ahat(m+j)=zhat(m+j)-app(phi,-1,zhat,m+j)+app(theta,-1,ahat,m+j);
X   end
X end
X k=len; % k is the number of unknown points
X for i=n:-1:1
X   if isnan(z(m+i)) 
X      lambda(m+i)=app(phi,1,lambda,m+i);
X      graduu(k)=ahat(m+i)+app(theta,1,lambda,m+i)-lambda(m+i);
X      k=k-1;
X   else
X      lambda(m+i)=ahat(m+i)+app(theta,1,lambda,m+i);
X   end
X end
X H=broyden(aest,aestvan,H,graduu,gradvan);
X gradvan=graduu;
X epsilon=graduu'*graduu;
end
X
SHAR_EOF
chmod 0644 inter2.m ||
echo 'restore of inter2.m failed'
Wc_c="`wc -c < 'inter2.m'`"
test 2748 -eq "$Wc_c" ||
	echo 'inter2.m: original size 2748, current size' "$Wc_c"
fi
# ============= inter2r.m ==============
if test -f 'inter2r.m' -a X"$1" != X"-c"; then
	echo 'x - skipping inter2r.m (File already exists)'
else
echo 'x - extracting inter2r.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'inter2r.m' &&
function [zhat,ahat]=inter2r(z,a,zest,phi,theta,tol)
X
% INTER2R finds estimates for missing observations in a time series.
% It is supposed that there are relatively few unknown observations,
% and that and ARIMA model for the series is already known.
% The call is of the form
%
%    [zhat,ahat]=inter2r(z,a,zest,phi,theta,tol)
%
% where
%   z     a time series with missing observations
%   a     estimates for the noise
%   zest  known values of the time series and initial estimates
%         for the missing observations
%   phi   the AR polynomial of the model
%   theta the MA polynomial of the model
%   tol  tolerance of the stopping criterion (squared norm of the gradient)
%   zhat  a time series including both the known values and 
%         estimates for the missing values
%   ahat  the noise estimates corresponding to zhat
X
% Method see the paper Ilkka Karanta: Optimization methods in time
% series interpolation. Communications in Statistics (to appear)
X
% Ilkka Karanta 1991-1993
X
% internal variables
%   lambda the co-state variable corresponding to the ARIMA equation
%   grad   the gradient of the cost function with respect to the
%          noise variables corresponding to missing observations
%   n      length of the part of time series with missing values
%   m      length of the part of time series with no missing values
X
apun=find(isnan(z));
m=apun(1)-1;
n=length(z)-m;
tyhjat=max(phi(length(phi),2), theta(length(theta),2));
zhat=z;
zhat(apun)=zest(apun);
ahat=a;
% the length of lambda is n+m+max(the order of phi, the order of theta)
lambda=zeros(n+m+tyhjat,1);
apu=z(find(isnan(z(m+1:m+n))));
H=eye(length(apu));
aest=a(apun);
graduu=zeros(length(apu),1);
gradvan=zeros(length(apu),1);
X
% initialization of the estimate vector
k=0;
pituus=length(apun);
for i=1:tyhjat
X  if ~isnan(z(i))
X    zhat(i)=z(i);
X  end
end
X
% initialization of ahat
for j=tyhjat+1:m+n
X   if isnan(z(j))
X      zhat(j)=zest(j);
X   else
X      zhat(j)=app(phi,-1,zhat,j)+ahat(j)-app(theta,-1,ahat,j);
X      ahat(j)=z(j)-zhat(j);
X      zhat(j)=z(j);
X   end
end
for j=m+n-tyhjat-1:-1:1
X   if isnan(z(j))
X      zhat(j)=zest(j);
X   else
X      zhat(j)=app(phi,1,zhat,j)+ahat(j)-app(theta,1,ahat,j);
X      ahat(j)=z(j)-zhat(j);
X      zhat(j)=z(j);
X   end
end
for j=tyhjat+1:m+n
X   if isnan(z(j))
X      zhat(j)=zest(j);
X   else
X      zhat(j)=app(phi,-1,zhat,j)+ahat(j)-app(theta,-1,ahat,j);
X      ahat(j)=z(j)-zhat(j);
X      zhat(j)=z(j);
X   end
end
% k=the number of missing observations
k=pituus;
for i=n:-1:1
X   if isnan(z(m+i))
X      lambda(m+i)=app(phi,1,lambda,m+i);
X      graduu(k)=ahat(m+i)+app(theta,1,lambda,m+i)-lambda(m+i);
X      k=k-1;
X   else
X      lambda(m+i)=ahat(m+i)+app(theta,1,lambda,m+i);
X   end;
end;
%for i=tyhjat+1:-1:1
%   lambda(i)=app(phi,1,lambda,i);
%   ahat(i)=lambda(i)-app(theta,1,lambda,i);
%end
sizeuu=size(graduu);
epsilon=graduu'*graduu;
X
% the actual loop
count=0;
pituus=length(apun);
X
while epsilon>tol
X count=count+1;
X aestvan=aest;
X aest=aest-H*graduu;
X ahat(1:tyhjat)=zeros(tyhjat,1);
X for i=1:length(apu)
X    ahat(apun(i))=aest(i);
X end
X for j=tyhjat+1:n+m
X   if isnan(z(j))
X      zhat(j)=app(phi,-1,zhat,j)+ahat(j)-app(theta,-1,ahat,j);
X   else
X      ahat(j)=zhat(j)-app(phi,-1,zhat,j)+app(theta,-1,ahat,j);
X   end
X end
X  for i=m+n-tyhjat-1:-1:1
X    if isnan(z(i))
X      zhat(i)=app(phi,1,zhat,i)+ahat(i)-app(theta,1,ahat,i);
X    else
X      ahat(i)=zhat(i)-app(phi,1,zhat,i)+app(theta,1,ahat,i);
X    end;
X  end
X for j=tyhjat+1:n+m
X   if isnan(z(j))
X      zhat(j)=app(phi,-1,zhat,j)+ahat(j)-app(theta,-1,ahat,j);
X   else
X      ahat(j)=zhat(j)-app(phi,-1,zhat,j)+app(theta,-1,ahat,j);
X   end
X end
X
X k=pituus; % k is the number of missing observations
X for i=n:-1:1
X   if isnan(z(m+i)) 
X      lambda(m+i)=app(phi,1,lambda,m+i);
X      graduu(k)=ahat(m+i)+app(theta,1,lambda,m+i)-lambda(m+i);
X      k=k-1;
X   else
X      lambda(m+i)=ahat(m+i)+app(theta,1,lambda,m+i);
X   end
X end
%for i=tyhjat+1:-1:1
%   lambda(i)=app(phi,1,lambda,i);
%   ahat(i)=lambda(i)-app(theta,1,lambda,i);
%end
X H=broyden(aest,aestvan,H,graduu,gradvan);
X gradvan=graduu;
X epsilon=graduu'*graduu;
end
X
SHAR_EOF
chmod 0644 inter2r.m ||
echo 'restore of inter2r.m failed'
Wc_c="`wc -c < 'inter2r.m'`"
test 4115 -eq "$Wc_c" ||
	echo 'inter2r.m: original size 4115, current size' "$Wc_c"
fi
# ============= intvar.m ==============
if test -f 'intvar.m' -a X"$1" != X"-c"; then
	echo 'x - skipping intvar.m (File already exists)'
else
echo 'x - extracting intvar.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'intvar.m' &&
function cova=intvar(phi,theta,sa2,z);
X
% INTVAR calculates the variances of errors of interpolated values of 
% missing observations. The estimates of missing observations are
% thought to be produced to minimize the least-squares criterion,
% and the time series is thought to be modelled by an ARIMA model.
% The call is of the form
%
% cova=intvar(phi,theta,sa2,z);
%
% cova  the covariance matrix of errors (the variances reside on the diagonal)
% phi   the AR polynomial of the time series
% theta the MA polynomial of the time series
% sa2   the variance of the residual noise
% z     the time series with missing values (missing values are
%       denoted by NaN)
X
% For a description of how the method works see the paper Ilkka Karanta:
% Optimization methods in time series interpolation (to appear in 
% Communications in Statistics)
X
% Ilkka Karanta 1991-1993
X
dimz=length(z);
apu=find(isnan(z));
% addbe is the number of missing values to be added before the 
%       beginning of the series to remove biasedness of the estimates
% adden the number of missing values to be added after the end of the 
%       series to remove biasedness of the estimates.
addbe=200;
adden=400;
inds=zeros(length(apu)+addbe+adden,1);
% additional unknown values are added before the beginning and after
% the end of the series to reduce the bias resulting from the effect
% of these values
inds(1:addbe)=1:addbe;
for i=1:length(apu) 
X  inds(addbe+i)=addbe+apu(i);
end
inds(length(apu)+addbe+1:length(apu)+addbe+adden)=addbe+dimz+1:addbe+...
X   dimz+adden;
dim=length(inds);
X
% inverse autocorrelations are calculated
[p,psi]=autocorr(theta,phi,dimz+addbe+adden+1);
X
ro(1:dimz+addbe+adden)=p(2:dimz+addbe+adden+1);
ro=ro / p(1);
X
Ri=zeros(dim);
for i=1:dim
X  for j=i+1:dim
X     Ri(i,j)=ro(inds(j)-inds(i));
X  end
end
X
Ri=Ri+Ri'+eye(dim);
X
cova=sa2/p(1)*inv(Ri);
X
X
SHAR_EOF
chmod 0644 intvar.m ||
echo 'restore of intvar.m failed'
Wc_c="`wc -c < 'intvar.m'`"
test 1858 -eq "$Wc_c" ||
	echo 'intvar.m: original size 1858, current size' "$Wc_c"
fi
# ============= pttopi.m ==============
if test -f 'pttopi.m' -a X"$1" != X"-c"; then
	echo 'x - skipping pttopi.m (File already exists)'
else
echo 'x - extracting pttopi.m (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'pttopi.m' &&
function pi=pttopi(phi,theta,n);
X
% PTTOPI calculates the ratio of two lag polynomials phi and theta.
% The call is of the form
%  pi=pttopi(phi,theta,n);
% where pi(B)=phi(B)/theta(B) (B is the lag operator)
%       phi is the numerator polynomial
%       theta is the denumerator polynomial
%       n     is the number of lags wanted in pi
%       Huom. length(pi)=n+1, koska ensimm{inen alkio pi(0) onkin pi(1)
X
% Reference: P.J. Brockwell and R.A. Davis: "Time series: 
% theory and methods" (Springer 1987), sec. 3.3
X
% Ilkka Karanta 1991-1993
X
phisize=size(phi);
thetasize=size(theta);
apuphi=zeros(n,1);
aputheta=zeros(n,1);
for i=1:phisize(1)
X  apuphi(phi(i,2))=phi(i,1);
end
for i=1:thetasize(1)
X  aputheta(theta(i,2))=theta(i,1);
end
X
for i=1:n
X  pi(i)=aputheta(i)-apuphi(i);
X  for j=1:i-1
X    pi(i)=pi(i)+aputheta(j)*pi(i-j);
X  end 
end
SHAR_EOF
chmod 0644 pttopi.m ||
echo 'restore of pttopi.m failed'
Wc_c="`wc -c < 'pttopi.m'`"
test 848 -eq "$Wc_c" ||
	echo 'pttopi.m: original size 848, current size' "$Wc_c"
fi
exit 0
