TIME SERIES MODELLING OF THE KOBE-OSAKA EARTHQUAKE RECORDINGS

A problem of great interest in monitoring a nuclear test ban treaty (NTBT) is related to interpreting properly the differences between a waveform generated by a nuclear explosion and that generated by an earthquake. With a view of comparing these two types of waveforms, Singh (1992) developed a technique for identifying a model in time domain. Fortunately this technique has been found useful in modelling the recordings of the killer earthquake occurred in the Kobe-Osaka region of Japan at 5.46 am on 17 January, 1995. The aim of the present study is to show how well the method for identifying a model (developed by Singh (1992)) can be used for describing the vibrations of the above mentioned earthquake recorded at Charters Towers in Queensland, Australia.


Introduction.
Many researchers have studied the structure of seismic records.Most of these studies (with the exception of Tjostheim [3,4] and Dargahi-Noubery et al. [2]) have been carried out in the frequency domain and they have been mainly concerned with the properties of their spectra.
The purpose of the present study is to model the recordings of the earthquake recently occurred at the Kobe-Osaka region of Japan on 17 January, 1995, in time domain, that is, to fit an autoregressive moving average model to the seismogram (see Figure 1.1) recorded at Charters Towers, Queensland, Australia.
One of the advantages of fitting a parametric model is that if a good model can be fitted, the process can be characterized by numerical values of a few parameters as opposed to the conjectural interpretation of the plots of spectra.Other advantages of the parametric modelling a time series in general have been highlighted by Box and Jenkins [1].
It has been found (see Tjostheim [3,4], Dargahi-Noubary et al. [2]) that the lowerorder autoregressive (AR) models are often appropriate for short period seismograms.Let Y (t) be a wide-sense-stationary (wss) process in discrete time.The process is said to be an AR process of order p (abbreviated AR(p)) if Y (t) satisfies the difference equation ϕ(B)Y (t) = ε(t), (1.1) where that E{Y (t)ε(s)} = 0, for s > t.Obviously, Y (t) has zero mean (an assumption valid for seismograms).Notice that model (1.1) contains p +1 parameters, namely, the ϕ i 's, i = 1,...,p and σ 2 ε .Once the order p is determined, the standard methods for estimating parameters can be used.
It is well known that a seismogram recorded over its entire range is zero-mean stationary in some sections while it is nonstationary-in-variance in others.
It may be noticed that the pattern of the series in Figure 1.1 is of this type.We divide the whole seismogram into three subseries and denote them by X 1 (t), X 2 (t), and X 3 (t) as shown in Figure 1.1.Thus we have three series out of the seismogram in Figure 1.1.Then a model of the following type may be fitted to the seismogram, namely (see Figure 1.1).Series X 1 (t), X 2 (t), and X 3 (t) consist of 670, 535, and 1051 observations, respectively, and each is assumed to follow an ARMA (p, 0) process.Now, the problem is to identify the model for each of these series and then to estimate the parameters.
Before discussing these problems in Sections 3 and 4, we review briefly the related work done previously in time domain, for information and ready reference to the reader.

Short review of previous studies of seismograms in time domain.
Tjostheim [3,4] considered the AR model defined by where X(t) denotes the seismic noise recorded at NORSAR, E(U(t)) = 0, Var(U(t)) = σ 2 u (t).It may be noted in model (2.1) that the coefficients α j (t), j = 1,...,p, and the residual variance σ 2 u (t) are assumed to be functions of t, though without any specific forms.They were estimated using two different samples and shown to be significantly different from the two samples indicating their time dependence.
However, in neither of these studies, nothing seems to have taken into account the fact that the series has increasing (or decreasing) variance in a certain range of the series.
It is this problem which has been addressed in this paper.A realistic approach is suggested in the following sections with an application to the Kobe-Osaka earthquake recordings.

3.
AR models with time-dependent error variance.Let (Ω,β,p) denote a probability space and let ᏸ(•) = ᏸ(Ω,β,p) denote the space of all real-valued random variables on (Ω,β,p) with zero mean and finite second-order moments.The space ᏸ(•) is called the Hilbert space if the inner product and norm are defined, respectively, by Consider now a zero-mean AR(P ) model with time-dependent error variance defined by where α(B) is the AR operator and w(t) ∈ ᏸ(•) is a white-noise with zero mean and time-dependent variance E{w where k and K are finite positive constants.Model (2.1) is stationary if, (i) all roots of equation α(x) = 0 lie outside the unit circle and (ii) the noise process Furthermore, we are interested in the family for which H w (t) is a parametric function of t.For all positive integer values of t, if the same parametric function is used, the class of possible functions will be restricted to the family such that H w (t) ≥ 0. Remark 3.1.It is interesting to note the relationship between Var X t = H X (t) and Var W t = H w (t) in the interval [t 1 ,t 2 ] as shown below: for convenience let t 1 = 0, t 2 = t, then for an AR(1), it can be seen that Notice that for |φ| < 1, and t → ∞, H X (t) will depend on the last two terms in (3.2) and hence the form of H X (t) will be similar to that of H w (t).For instance if H w (t) is a linear function of t, H X (t) will also be a linear function of t, approximately.
Remark 3.3.It may be noticed from Theorem 3.2(iii) that the correlation structure of X(t) is the same as that of the zero-mean stationary process {Y (t)}; although X(t) is zero-mean nonstationary-in-variance.
It is this common characteristic of the two processes which has been exploited in the sequel for the identification of both X(t) and Y (t) as well as the identification of the covariance structure of the white-noise w(t) in (3.2).
For example, let {Y (t)} be a zero-mean stationary AR(2) process defined by Under the transformation X(t) = S(t)Y (t), we have where In general, it can be seen that lim In particular, let S(t) = exp(α t), then we have  1),...,X(200)} was generated using (3.16).These samples along with their respective ACF's and PACF's are plotted in Figures 3.1a and 3.1b.

Inference and conclusions.
It may be noticed from Figure 3.1 that (i) the series Y (t) is stationary as expected; (ii) the series X(t) under the transformation (3.16) is a naturally exponentially increasing series since the variance of the white-noise w(t) is assumed to be an exponential function of t; (iii) the ACF's and PACF's of both Y (t) and X(t) are the same as expected (see Theorem 3.2), except for some sampling fluctuations; (iv) both ACF and PACF of X(t) (and/or Y (t)) suggest that the underlying stationary process Y (t) is an AR( 2) and the oscillatory nature of the ACF of either X(t) or Y (t), that one of the coefficients of the process Y (t), is negative.

Preliminary identification.
Given a zero-mean and nonstationary-in-variance series such as X(t) in the preceding section, the following steps are suggested to identify the process: (1) plot the series, its ACF and PACF, (2) if the plot of the series is exponentially increasing (or decreasing) on both sides of the mean then it may be argued that the error-variance is an exponential function of time t with positive (or negative) exponent.Similar interpretation may be given to the series if it is linearly increasing (or decreasing), (3) the ACF and PACF together would suggest the order of the underlying stationary process such as the Y (t) in the above example, (4) thus having identified the order of the AR process that can be fitted to the given X(t) series and the covariance structure of the white-noise, any standard procedure of estimation such as the maximum likelihood (ML) can be used to estimate the parameters in the model and the covariance function of the whitenoise.

Advantages.
One of the advantages of the above procedure is that it enables us to estimate the variance of the white-noise and the coefficients of the underlying stationary process which is referred to as Y (t) in the preceding example.For instance, let Y (t) be the underlying stationary AR(1) defined by and let then where Assuming S(t) = exp(αt), we have where β = σ 2 ε and c = 2α.Given X(1),...,X(n), the ML method will yield the estimators of ϕ * , β and thus enabling us to estimate σ 2 ε , the coefficients of models Y (t) and X(t) in (4.1) and (4.3), respectively.

5.
Fitting models to the Kobe-Osaka seismic recordings.The seismogram of the Kobe earthquake recorded at Charters Towers is given in Figure 1.1.The graph displays the vertical component of the ground displacement (in nm) versus time in minutes.The total number of observations recorded were 2256.This entire set of observations has been divided into three main parts.We denote them by X 1 (t), X 2 (t), and X 3 (t) series as shown in Figure 1.1.These series consist of 670, 535, and 1051 observations, respectively.

Modelling of series X 2 (t).
First of all, we study and model series X 2 (t) due to its unusual appearance.This consists of 535 values and it is displayed in Figure 5.1.
The ACF and PACF of X 2 (t) are shown in Figure 5.2.
Inference.(i) From Figure 5.3, we infer that the variance structure of the series X 2 (t) is time-dependent and seems to be an exponential function of t multiplied by a constant which may be taken as an error-variance of the underlying process Y 2 (t).(See Section 4.) (ii) Both ACF and PACF suggest that the representative process is pure autoregressive process of order three (R(3)).
Based on this preliminary inference, we postulate the underlying stationary process of the type where ε 3 (t) ∼ N(0,σ 2 ε ), and select S(t) = exp(αt).Put (5.2) Then the model for X 2 (t) may be defined by where then where ( Thus based on (i) and (ii) above we may infer that the modelling of series X 2 (t), following the procedures outlined in Sections 3 and 4, is more than satisfactory.
In the following subsection, we consider the modelling series X 3 (t).From Figures 5.5 and 5.6, it appears that the series in question represents a pure AR(4) process.Consequently, using ITSM, we fitted three models to X 3 (t), namely, AR(3), AR(4), and AR (5) and obtained the respective AICC-values along with the corresponding maximum likelihood estimates.The AICC-values are shown in Table 5.1.
The series of X 3 (t) consisting of 1051 observations is plotted in Figure5.5.Its ACF and PACF are plotted inFigures 5.6a and 5.6b, respectively.
(a) Plot of the ACF of X2 (t).

Figure 5 . 4 X 3 (t) 0 Figure 5 . 5
Figure 5.4 Figure 5.6 On comparing the plots of the original series X2 (t), its ACF, and PACF in Figures 5.1 and 5.2 with the corresponding plots of the estimated series X2 (t), its ACF and PACF in Figures 5.3, 5.4a, and 5.4b, one may note the striking similarity.