A new family of time series models, called the Full Range Autoregressive model, is introduced which avoids the difficult problem of order determination in time series analysis. Some of the basic statistical properties of the new model are studied. Further, the paper describes the Bayesian inference and forecasting as applied to the Full Range Autoregressive model. The Canadian lynx data is used to compare the efficiency of the predictive power of the new model with those of some of the existing models in the time series literature.

Keywords: full range autoregressive model, identifiability, stationary condition, posterior distribution, bayesian predictive distribution

 The early days of time series analysis, most of the models fitted to the real life data were restricted to low orders because of availability of high speed computers and other facilities. However, now with the availability of high speed computers, there is no need for this type of restriction on the order determination and estimation of the fitted models. Further, most of the work in time series analysis are concerned with series having the property that the degree of dependence between observations, separated by a long time span, is zero or highly negligible. However, the empirical studies by Lawrance and Kottegoda¹ reveal, particularly in cases arising in economics and hydrology, that the degree of dependence between observations a long time span apart, though small, is by no means negligible. Therefore, there is still a need for a family of models which can fully depict the properties of stationarity, linearity and long range dependence.

Moreover, the existing theory of autoregressive models assume that the coefficients of the model are not connected in any way among each other. Therefore, it would be useful, from practical point of view, to propose new models, called the Full Range Auto Regressive model and denoted as FRAR model for short, which can accommodate long range dependence and have the property that the coefficients of the past values in the model are functions of a limited number of parameters.

Thus, the chief objective of this paper is to introduce a family of new models which would involve only a few parameters and at the same time incorporate long range dependence, which would be an acceptable alternative to the current models representing stationary time series.

A family of models, introduced in this paper, called Full Range Auto Regressive model and denoted as FRAR model for short, are defined in such a way that they possess the following basic features. 

1. The models should be capable of representing long term persistence. This is justified by the fact that the future may not depend on the present and a few past values alone, but may depend on the present and the whole past.
2. The parameters of the model, which are likely to be large in number due to (1), should exhibit some degree of dependence among themselves.

Therefore, the new models are expected to have infinite structure with a finite number of parameters and so completely avoid the problem of order determination.

An outline of this paper is as follows. In Section 2, the FRAR model is defined, the identifiability region is obtained, the stationarity condition is derived, and the asymptotic stationarity is studied. In Section 3, the Bayesian analysis of the FRAR model is discussed and the predictive density of a single future observation is derived. In Section 4 the Canadian lynx data is used for forecasting through the FRAR model. In Section 5 a comparative study is provided to examine the efficiency of FRAR model. In Section 6 the summary and conclusion is given.

The model

We define a family of models by a discrete-time stochastic process $\left(X_t\right)$ , $t=0,\pm 1,\pm 2,\mathrm{...}$ , called the Full Range Auto Regressive (FRAR) model, by the difference equation

$X_t={\displaystyle {\sum }_{r=1}^{\infty }{a}_r{X}_{t-r}+e_t}$       (1)

where $a_r=k\sin \left(r\theta \right)\cos \left(r\varphi \right)/{\alpha }^r$ , $\left(r=1,2,\text{ 3,}\mathrm{...}\right)$ , $k$ , $\alpha$ , $\theta$  and $\varphi$  are parameters, $e_1$ , $e_2$ , $e_3$ , … are independent and identically distributed normal random variables with mean zero and variance ${\sigma }^2$ . The initial assumptions about the parameters are as follows:

It is assumed that $X_t$  will influence $X_{t+n}$  for all positive n and the influence of $X_t$  on $X_{t+n}$  will decrease, at least for large n, and become insignificant as n becomes very large, because more important for the recent observations and less important for an older observations. Hence $a_n$  must tend to zero as n goes to infinity. This is achieved by assuming that $\alpha >1$ . The feasibility of $X_t$  having various magnitudes of influence on $X_{t+n}$ , when n is small, is made possible by allowing k to take any real value. Because of the periodicity of the circular functions sine and cosine, the domain of $\theta$  and $\varphi$  are restricted to the interval $\left[0,2\pi \right)$ .

Thus, the initial assumptions are $\alpha >1$ , $k\in R$ , and $\theta$ , $\varphi \in \left[0,2\pi \right)$ . i.e., $\Theta =\left(\alpha ,k,\theta ,\varphi \right)\in S*$ , where $S*=\left\{\alpha ,k,\theta ,\varphi |k\in R,\alpha >1,\theta ,\varphi \in \left[0,2\pi \right)\right\}$ . Further restrictions on the range of the parameters are placed by examining the identifiability of the model.

Identifiability condition

Identifiability ensures that there is a one to one correspondence between the parameter space and set of associated probability models. Without identifiability it is meaningless to proceed to estimate the parameters of a model using a set of given data. In the present context, identifiability is achieved by restricting the parameters space in such a way that no two points in the parameter space could produce the same time series model.

The coefficients $a_n$ ’s in (1) are functions of $k$ , $\alpha$ , $\theta$ , $\varphi$  as well as n. That is, $a_n=a_n\left(k,\alpha ,\theta ,\varphi \right)=k\sin \left(n\theta \right)\cos \left(n\varphi \right)/{\alpha }^n$ , $\theta \in S*$ , $n=1,2,\text{ 3,}\mathrm{...}$ .

Define $A=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\pi \le \theta ,\varphi <2\pi \right\}$ ,

$B=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\text{ 0}\le \theta <\pi ,\pi \le \varphi <2\pi \right\}$ , (2)

$C=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\pi \le \theta <2\pi ,0\le \varphi <\pi \right\}$ ,

$D=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,0\le \theta ,\varphi <\pi \right\}$ .

Since $a_n=a_n\left(k,\alpha ,\theta ,\varphi \right)=a_n\left(-k,\alpha ,\text{ 2}\pi -\theta ,\text{ 2}\pi -\varphi \right)$ , $\theta \in S*$           

to each $\left(\alpha ,k,\theta ,\varphi \right)$  belonging to A there is $a$   $\left(\alpha ,k,\theta \text{'},\varphi \text{'}\right)\left(\theta \text{'}=2\pi -\theta \text{ and }\varphi \text{'}=2\pi -\varphi \right)$  belonging to D such that $a_n\left(k,\alpha ,\theta ,\varphi \right)=a_n\left(-k,\alpha ,\theta \text{'},\varphi \text{'}\right)$ . So A is omitted. Similarly, it can be shown that B and C can also be omitted.
Define    $D_1=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\pi \text{/2}\le \theta ,\varphi <\pi \right\}$ ,

$D_2=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\text{ 0}\le \theta <\pi \text{/2},\pi \text{/2}\le \varphi <\pi \right\}$ ,

$D_3=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\text{ 0}\le \theta ,\varphi <\pi \text{/2}\right\}$ ,

$D_4=\left\{\alpha ,k,\theta ,\varphi |\alpha >1,k\in R,\pi \text{/2}\le \theta <\pi ,\text{ 0}\le \varphi <\pi \text{/2}\right\}$ .

Since $a_n\left(k,\alpha ,\theta ,\varphi \right)=a_n\left(-k,\alpha ,\pi -\theta ,\pi -\varphi \right)$

 for $k\in R,\alpha >1,\text{ 0}\le \theta ,\varphi <\pi$                                                (3)

Using (3) it can be shown as before, that the regions $D_1$  and $D_2$  can be omitted. Since no further reduction is possible, it is finally deduced that the region of identifiability of the model is given by $S=\left\{\alpha ,k,\theta ,\varphi |k\in R,\alpha >1,\theta \in \left[0,\pi \right),\varphi \in \left[0,\pi /2\right)\right\}$ .

Stationarity of the FRAR process

The stationarity of the newly developed FRAR time series model is now examined. The model is given by $X_t={\displaystyle {\sum }_{r=1}^{\infty }{a}_r{X}_{t-r}+e_t}$ . That is, $\left(1-a_{1}B-a_2{B}^2-\dots \right)X_t=e_t$ , where B is the backward shift operator, defined by $B^n{X}_t=X_{t-n}$ . Thus, the model is given by $\Psi \left(B\right)X_t=e_t$ , or $X_t={\Psi }^{-1}\left(B\right)e_t$ , where $\Psi \left(B\right)=1-a_{1}B-a_2{B}^2-\dots$ .

Box and Jenkins² and Priestley³ have shown that a necessary condition for the stationarity of such processes is that the roots of the equation $\Psi \left(B\right)=0$  must all lie outside the unit circle. So, it is now proposed to investigate the nature of the zeros of $\Psi \left(B\right)$ .

The power series $\Psi \left(B\right)$  may be rewritten as $\Psi \left(B\right)=1-\left[a_{1}B+a_2{B}^2+\dots \right]=1-{\displaystyle {\sum }_{n=1}^{\infty }{a}_n{B}^n}$ , where $a_n{B}^n=\left(k{B}^n/{\alpha }^n\right)\left[\sin \left(n\theta \right)\cos \left(n\varphi \right)\right]=\left(k\text{'}{B}^n/{\alpha }^n\right)\left[\sin \left(n{\theta }_1\right)+\sin \left(n{\theta }_2\right)\right]$ , $k\text{'}=k/2$ , ${\theta }_1=\theta +\varphi$  and ${\theta }_2=\theta -\varphi$ . Therefore, ${\displaystyle {\sum }_{n=1}^{\infty }{a}_n{B}^n}={\displaystyle {\sum }_{n=1}^{\infty }\frac{k\text{'}{B}^n}{{\alpha }^n}\sin \left(n{\theta }_1\right)+}{\displaystyle {\sum }_{n=1}^{\infty }\frac{k\text{'}{B}^n}{{\alpha }^n}\sin \left(n{\theta }_2\right)}$ .

The above two series are separately evaluated below.

${\displaystyle {\sum }_{n=1}^{\infty }\frac{k\text{'}{B}^n}{{\alpha }^n}\sin \left(n{\theta }_1\right)=IP\text{ of }}{\displaystyle {\sum }_{n=1}^{\infty }\frac{k\text{'}{B}^n}{{\alpha }^n}{e}^{in{\theta }_1}=IP\left\{kB{e}^{i{\theta }_1}{\left(\alpha -B{e}^{i{\theta }_1}\right)}^{-1}\right\}}=k\text{'}B\alpha \sin \left({\theta }_1\right)/G_1$ , where $G_1=B^2+{\alpha }^2-2B\alpha \cos \left({\theta }_1\right)$ and IP stands for imaginary part.

Similarly, it can be shown that ${\displaystyle {\sum }_{n=1}^{\infty }\frac{k\text{'}{B}^n}{{\alpha }^n}\sin \left(n{\theta }_2\right)=\left(k\text{'}B\alpha \sin \left(n{\theta }_2\right)\right)/G_2}$ , where $G_2=B^2+{\alpha }^2-2B\alpha \cos \left({\theta }_2\right)$ .

Therefore,

${\displaystyle {\sum }_{n=1}^{\infty }{a}_n{B}^n}=k\text{'}B\alpha \left[\left(B^2+{\alpha }^2\right)\left(\sin \left({\theta }_1\right)+\sin \left({\theta }_2\right)\right)-2B\alpha \left(\sin \left({\theta }_1\right)\cos \left({\theta }_2\right)-\cos \left({\theta }_1\right)\sin \left({\theta }_2\right)\right)\right]/G_1{G}_2$ .

Thus, $\Psi \left(B\right)=1-{\displaystyle {\sum }_{n=1}^{\infty }{a}_n{B}^n}=0$ implies that $\left(B^2+{\alpha }^2-2B\alpha \cos \left({\theta }_1\right)\right)\left(B^2+{\alpha }^2-2B\alpha \cos \left({\theta }_2\right)\right)-k\text{'}B\alpha {\left[\left(B^2+{\alpha }^2\right)s_1-2B{d}_1{c}_2\right]}^{}=0$

where $c_1=\cos \left({\theta }_1\right)+\cos \left({\theta }_2\right)=2\cos \left(\theta \right)\cos \left(\phi \right)$ , $c_2=\sin \left(2{\theta }_2\right)$ , $s_1=\sin \left({\theta }_1\right)+\sin \left({\theta }_2\right)=2\sin \left(\theta \right)\cos \left(\phi \right)$ , $d_1=\cos \left({\theta }_1\right)-\cos \left({\theta }_2\right)$ . After simplifying, the above equation becomes $B^4-B^3\alpha \left(2c_1+k\text{'}{s}_1\right)+B^2{\alpha }^2\left(2+4d_1+2k\text{'}{c}_2\right)-B{\alpha }^3\left(2c_1+k\text{'}{s}_1\right)+{\alpha }^4=0$ . Thus,

$B^4-B^3\alpha A_1+B^2{\alpha }^2{A}_2-B{\alpha }^3{A}_1+{\alpha }^4=0$ (4)

or $S^4-A_1{S}^3+A_2{S}^2-A_{1}S+1=0$ (5)

where $A_1=2c_1+k\text{'}{s}_1=\cos \left(\phi \right)\left(4\cos \left(\theta \right)+k\sin \left(\phi \right)\right)$ , $A_2=2+4d_1+2k\text{'}{c}_2=2\left[1-\sin \left(\varphi \right)\left(4\sin \left(\theta \right)-k\cos \left(\varphi \right)\right)\right]$ , and $S=B/\alpha$ . This equation (of degree 4) reduces to $Z^2-A_{1}Z+\left(A_2-2\right)=0$  where $Z=S+\left(1/S\right)$ .

The roots of this equation are, say $r_1$  and $r_2$ , are given by $Z=\left(1/2\right){\left[A_1\pm \sqrt{\left(A_1^2-4A_2+8\right)}\right]}^{}$ .

Since $Z=S+\left(1/S\right)$ , one finally gets the four roots of the equation (4), as $R_1=\left(1/2\right){\left[r_1+\sqrt{r_1^2-4}\right]}^{}$ , $R_2=\left(1/2\right){\left[r_1-\sqrt{r_1^2-4}\right]}^{}$ , $R_3=\left(1/2\right){\left[r_2+\sqrt{r_2^2-4}\right]}^{}$  and $R_2=\left(1/2\right){\left[r_2-\sqrt{r_2^2-4}\right]}^{}$ .

The equation (5) implies that, if $S_0$  is a root of the equation (5) then $1/S_0$  is also a root. This implies that $\alpha S_0$  and $\left(\alpha /S_0\right)$  are roots of equation (4). Therefore the process is stationary for sufficiently large values of $\alpha$ . But when $\alpha$  is small it seems difficult to examine the stationarity of the process by this approach. Hence, it is proposed to study the asymptotic stationarity of the process in the following section.

Asymptotic stationarity of the FRAR process

In this section we derive the condition for asymptotic stationarity of the FRAR process. For which one has to solve the difference equation (1), so as to obtain an expression for $X_t$  in terms of $e_t$ , $e_{t-1}$ , $e_{t-2}$ , $e_{t-3}$ , .... The precise solution of this equation depends on the initial conditions. So to investigate the nature of the first and second moments of $X_t$ , following Priestley,³ it is assumed that $X_t=0$  for $t<-N$ , N being the number of observations in the time series. Then solving (1) by repeated substitutions one obtains

$X_t=e_t+a_{11}{X}_{t-1}+a_{12}{X}_{t-2}+a_{13}{X}_{t-3}+\dots$ ,

where $a_{1j}=a_j=\left(k/{\alpha }^j\right)\sin \left(j\theta \right)\cos \left(j\varphi \right)$ ; $j=1,2,\dots$ ,

$=e_t+a_{11}{e}_{t-1}+a_{22}{X}_{t-2}+a_{23}{X}_{t-3}+a_{24}{X}_{t-4}\dots$ ,

where $a_{2j}=a_{11}{a}_{1j-1}+a_{1j}$ ; $j=2,3,4\dots$ .

Similarly proceeding one finally gets

$X_t={\left[e_t+a_{11}{e}_{t-1}+a_{22}{e}_{t-2}+a_{33}{e}_{t-3}+a_{44}{e}_{t-4}+\dots +a_{pp}{e}_{t-p}\right]}^{}$ $+{\left[a_{p+1p+1}{X}_{t-\left(p+1\right)}+a_{p+1p+2}{X}_{t-\left(p+2\right)}+\dots \right]}^{}$  

where $a_{ij}=a_{i-1i-1}{a}_{1j+1-i}+a_{i-1j}$  with $j>i=2,3,4\dots$ . Thus, if it is assumed that $X_t=0$  for $t\le -N$ , which implies has $n=N+t-1$ , then, $X_t=e_t+a_{11}{e}_{t-1}+a_{22}{e}_{t-2}+a_{33}{e}_{t-3}+a_{44}{e}_{t-4}+\dots +a_{N+t-1,N+t-1}{e}_{1-N}$ .

Further, it can be shown that

$E\left[X_t{X}_{t+1}\right]={\sigma }_e^2{\left[a_{11}\left(1+a_{11}^2+a_{22}^2+\dots +a_{N+t-1N+t-1}^2\right)+\left(a_{11}{a}_{12}+a_{22}{a}_{23}+\dots +a_{N+t-2N+t-2}{a}_{N+t-2N+t-1}\right)\right]}^{}$

$E\left[X_t{X}_{t+2}\right]={\sigma }_e^2{[a_{22}\left(1+a_{11}^2+a_{22}^2+\dots +a_{N+t-1N+t-1}^2\right)}^{}$

$+a_{11}\left(a_{11}{a}_{12}+a_{22}{a}_{23}+\dots +a_{N+t-2N+t-2}{a}_{N+t-2N+t-1}\right)$

$+\left(a_{11}{a}_{13}+a_{22}{a}_{24}+\dots +a_{N+t-3N+t-3}{a}_{N+t-3N+t-1}\right)]$

$E\left[X_t{X}_{t+3}\right]={\sigma }_e^2{[a_{33}\left(1+a_{11}^2\left(1+a_{22}+a_{33}+\dots +a_{N+t-3N+t-3}^2\right)\right)}^{}$

$+a_{11}\left(a_{22}{a}_{34}+a_{33}{a}_{45}+\dots +a_{N+t-3N+t-3}{a}_{N+tN+t-1}\right)$

${+\left(a_{11}{a}_{34}+a_{22}{a}_{45}+\dots +a_{N+t-1N+t-1}{a}_{N+t-3N+t-3}\right)]}^{}$  

and in general

$E\left[X_t{X}_{t+s}\right]={\sigma }_e^2\left[a_{ss}+a_{11}{a}_{s+1s+1}+a_{22}{a}_{s+2s+2}+\dots +a_{N+t-1N+t-1}{a}_{N+t+s-1N+t+s-1}\right]$  

where $a_{ss}=a_{11}{a}_{s-1s-1}+a_{s-1s}$ . Therefore, allowing $N\to \infty$ , we get $E\left[X_t\right]=0$ , $Var\left[X_t\right]={\sigma }_e^2\left[1+a_{11}^2+a_{22}^2+\dots \right]$  and $E\left[X_t{X}_{t+s}\right]={\sigma }_e^2{\left[a_{ss}+a_{11}{a}_{s+1s+1}+\dots \right]}^{}$  provided the series on the right converges. Thus, it is seen that if $E\left[X_t{X}_{t+s}\right]$  exists then it is a function of s only. In order to examine the convergence of $Var\left[X_t\right]$  and $E\left[X_t{X}_{t+s}\right]$ , first the behaviour of $a_{ij}$ , as j tends infinity, is investigated. Since $a_{1j}=a_j=\left(k/{\alpha }^j\right)\sin \left(j\theta \right)\cos \left(j\phi \right)$ , $\left|a_{1j}\right|\le \left|k\right|/{\alpha }^j$ . Similarly, $\left|a_{2j}\right|\le \left|k\right|\left(1+\left|k\right|\right)/{\alpha }^j$ ; $j\ge 2$ . Thus, in general $\left|a_{nj}\right|\le \left|k\right|{\left(1+\left|k\right|\right)}^{n-1}/{\alpha }^j$ , for $j\ge n$ .

Since $\alpha >1$ , the above relation implies that $\left|a_{nj}\right|\to 0$  as $j\to \infty$ , for any fixed n. Thus ${\displaystyle {\sum }_{n=1}^{\infty }{a}_{jj}^2}$  will converge if $\left|\frac{\left(1+\left|k\right|\right)}{\alpha }\right|<1$ .

If we assume that $1-\alpha <k<\alpha -1$ , then one can show that $Var\left[X_t\right]={\sigma }_{X_t}^2\le {\sigma }_e^2\frac{k^2}{{\left(1+k\right)}^2}\left[\frac{{\alpha }^2}{{\alpha }^2-{\left(1+k\right)}^2}\right]$  and $E\left[X_t{X}_{t+s}\right]\le {\sigma }_e^2\frac{k^2}{{\left(1+k\right)}^2}\frac{k{\left(1+k\right)}^{s-1}}{{\alpha }^s}\left[\frac{{\alpha }^2}{{\alpha }^2-{\left(1+k\right)}^2}\right]$ .

Therefore, the auto-correlation function of the process exists and, as shown earlier, it is a function of s only. Finally allowing $t\to \infty$ , it is seen that

1. $\underset{t\to \infty }{\lim }E\left[X_t\right]$  and $\underset{t\to \infty }{\lim }Var\left[X_t\right]$  exist finitely;
2. $\underset{t\to \infty }{\lim }Cov\left[X_t,X_{t+s}\right]$  exists finitely and is a function of ‘s’ only.

Thus, the condition for $\left\{X_t\right\}$  to be asymptotically stationary is that $1-\alpha <k<\alpha -1$ . Therefore, we summarized the above results by the following theorem 1.

Theorem 1: The Full Range Auto Regressive (FRAR) process $\left\{X_t\right\}$  is asymptotically stationary and identifiable if and only if the domain of the parameter space $S$ is

$\left\{k\text{, }\alpha \text{, }\theta \text{, }\varphi /k\in R,1-\alpha <k<\alpha -1,\theta \in \left[0,\pi \right),\varphi \in \left[0,\pi /2\right)\right\}$ $\alpha >1$ .

Thus, the new FRAR model incorporates long range dependence, involves only four parameters and is totally free from order determination problems.

The posterior analysis

The Bayesian approach to the analyses of the new model consists in determining the posterior distribution of the parameters of the FRAR model and the predictive distribution of future observations. From the former, one makes posterior inferences about the parameters of the FRAR model including the variance of the white noise. From the latter, one may forecast future observations. All these techniques are illustrated by Broemeling⁴ for autoregressive models.

We shall consider the FRAR model and assume that it is asymptotically stationary and identifiable.

The problem is to estimate the unknown parameters $k$ , $\alpha$ , $\theta$ , $\varphi$  and ${\sigma }^2$ , using the Bayesian methodology on the basis of a past random realization of $\left\{X_t\right\}$  say $x=\left(x_1,x_2,\mathrm{...},x_N\right)$ .

The joint probability density of $X_1,X_2,\mathrm{...},X_N$  is given by

$P\left(X/\Theta \right)\propto {\left({\sigma }^2\right)}^{-N/2}\exp \left[-\frac{1}{2{\sigma }^2}{\displaystyle {\sum }_{t=1}^2{\left(x_t-k{\displaystyle {\sum }_{r=1}^{\infty }{a}_r{x}_{t-r}}\right)}^2}\right]$  (6)

where $x=\left(x_1,x_2,\mathrm{...},x_N\right)$ , $\Theta =\left(k,\alpha ,\theta ,\phi ,{\sigma }^2\right)$ and $a_r=\left(1/{\alpha }^2\right)\sin \left(r\theta \right)\cos \left(r\varphi \right)$ .

The notation $P$  is used as a general notation for the probability density function of the random variables given within the parentheses following $P$ and $X_0,X_{-1},X_{-2},\mathrm{...}$  are the past realizations on $X_t$  which are unknown. Following Priestley² and Broemeling,³ these are assumed to be zero for the purpose of deriving the posterior distribution of $\Theta$ . Therefore, the range for the index r, viz., 1 through ∞, reduces to 1 through N and so, in the joint probability density function of the observations given by (6), the range of the summation 1 through ∞ can be replaced by 1 through N. By expanding the square in the exponent and simplifying, one gets
$P\left(X/\Theta \right)\propto {\left({\sigma }^2\right)}^{-N/2}\exp (-Q/2{\sigma }^2)$ (7)

where $Q=T_{00}+k^2{\displaystyle {\sum }_{r=1}^N{a}_r^2{T}_{rr}}+2k^2{\displaystyle {\sum }_{r<s;r,s=1}^N{a}_r{a}_s{T}_{rs}}-2k{\displaystyle {\sum }_{r=1}^N{a}_r{T}_{r0}}$ , $T_{rs}={\displaystyle {\sum }_{t=1}^N{x}_{t-r}{x}_{t-s}}$ , $r,s=0,1,\mathrm{...},N$ , $\Theta \in S$ .

To find the posterior distribution of $\Theta$ we first have to specify the prior distribution for the parameters.

$\alpha$ is distributed as the displaced exponential distribution(since it is bigger than 1) with parameter $\beta$ ; ${\sigma }^2$  has the inverted gamma distribution with parameter v and δ; $k$ , $\theta$  and $\varphi$  are uniformly distributed over their domain.

Thus, the joint prior density function of $\Theta$   $\left(\Theta \in S\right)$ is given by

$P\left(\Theta \right)\propto \beta \exp \left(-\beta \left(\alpha -1\right)-\nu /{\sigma }^2\right){\left({\sigma }^2\right)}^{-\left(\delta +1\right)}$          (8)

$P\left(\Theta /X\right)\propto {\left({\sigma }^2\right)}^{-N/2}\exp (-Q/2{\sigma }^2)\exp {\left[-\beta (\alpha -1)-\nu /{\sigma }^2\right]}^{}{({\sigma }^2)}^{-(\delta +1)}$                     (9)

$\propto \exp \left[-\beta (\alpha -1)\right]\exp {\left[-1/2{\sigma }^2(Q+2\nu )\right]}^{}{({\sigma }^2)}^{-\left[\left(N/2\right)+\delta +1\right]}$                     (10)

Integrating ${\sigma }^2$ out of this joint posterior distribution, we obtain the joint posterior distribution of $k$ , $\alpha$ , $\theta$  and $\varphi$ ,

$P\left(k,\alpha ,\theta ,\varphi /X\right)\propto e^{-\beta (\alpha -1)}{\left\{C\left[1+A_1{\left(k-B_1\right)}^2\right]\right\}}^{-d}$                                                                 (11)

Thus, the posterior distribution of k conditional on $\alpha$ , $\theta$  and $\varphi$  is a t-distribution located at $B_1$  with $\left(2d-1\right)$  degrees of freedom.

Thus, the joint posterior density function of $\alpha$ , $\theta$  and $\varphi$  can be obtained by integrating with respect to k. Thus,

$P\left(\alpha ,\theta ,\varphi /x\right)\propto \exp \left(-\beta (\alpha -1)\right)C^{-d}{A}_1^{-1/2}$ ; with $\alpha >1$ , $0\le \theta <\pi$  and $0\le \varphi <\pi /2$ .      (12)

The above joint posterior density of $\alpha$ , $\theta$  and $\varphi$  is a very complicated expression and is analytically intractable. One way of solving the problem is to find the marginal posterior density of $\alpha$ , $\theta$  and $\varphi$  from the joint density (12) using ordinary numerical integration, using FORTRAN.

One-step-ahead prediction

In order to forecast $x_{N+1}$  using the random realization $x_1,x_2,\mathrm{...},x_N$  on $\left(X_1,X_2,\mathrm{...},X_N\right)$ , one must find the conditional distribution of $X_{N+1}$  given the past observations. This is the predictive distribution of $X_{N+1}$  and will be derived by multiplying the conditional density of $X_{N+1}$  given $X_1,X_2,\mathrm{...},X_N$ , $\Theta$  and the posterior density of $\Theta$  given $X_1,X_2,\mathrm{...},X_N$  and then integrating with respect to $\Theta$ . That is, $P\left(X_{N+1}/X_1,X_2,\mathrm{...},X_N\right)={\displaystyle {\int }_{\Theta }P\left(X_{N+1}/X_1,X_2,\mathrm{...},X_N,\Theta \right)P\left(\Theta /X_1,X_2,\mathrm{...},X_N\right)d\Theta }$ .

Thus, we obtain

$P\left(x_{N+1}/x_1,x_2,\mathrm{...},x_N,\Theta \right)\propto {\left({\sigma }^2\right)}^{-1/2}\exp \left[-\frac{1}{2{\sigma }^2}{\left(x_{N+1}-k{\displaystyle {\sum }_{i=1}^{\infty }{a}_i{x}_{N+1-i}}\right)}^2\right]$ , $x_{N+1}\in R$ . (13)

The square in the exponent in the above expression, say $Q_1$ , can be rewritten, after expanding the square, as $Q_1=x_{N+1}^2+k^2{\displaystyle {\sum }_{i=1}^N{a}_i^2{P}_i^2+}2k^2{\displaystyle {\sum }_{i<j;i=1}^N{a}_i^{}{a}_j^{}{P}_{ij}^{}-2k}{\displaystyle {\sum }_{i=1}^N{a}_i^{}{P}_i^{}{X}_{N+1}}$ , where $P_i=X_{N+1-i}$  and $P_{ij}=X_{N+1-i}{X}_{N+1-j}$ . Now multiplying (13) by the joint posterior density of $\Theta$  and integrating over the parameter space $\Theta$ , we obtain,

$P\left(x_{N+1}/x_1,x_2,\mathrm{...},x_N\right)\propto {\displaystyle \int \exp \left(-\beta \left(\alpha -1\right)\right){\left(1/{\sigma }^2\right)}^{\left[\frac{N}{2}+\delta +\frac{1}{2}+1\right]}}\exp \left[-\frac{1}{2{\sigma }^2}\left(Q+Q_1+2\upsilon \right)\right]d\Theta$ (14)

First, integrating out ${\sigma }^2$  in (14), one gets the joint distribution of $x_{N+1}$ , $k$ , $\alpha$ , $\theta$  and $\varphi$  as

$P\left(x_{N+1},k,\alpha ,\theta ,\varphi /x_1,x_2,\mathrm{...},x_N\right)\propto \exp \left(-\beta \left(\alpha -1\right)\right){\left(Q+Q_1+2\upsilon \right)}^{-\left(\frac{N+1}{2}+\delta \right)}$                            (15)

where $d_1={\displaystyle {\sum }_{i=1}^N{a}_i^2{T}_{ii}+}2{\displaystyle {\sum }_{i<j;i=1}^N{a}_i{a}_j{T}_{ij}}$ , $d_2={\displaystyle {\sum }_{i=1}^N{a}_i^2{P}_i^2+}2{\displaystyle {\sum }_{i<j;i=1}^N{a}_i{a}_j{P}_{ij}}$ , $d_3={\displaystyle {\sum }_{i=1}^N{a}_i{T}_{i0}}$ , $d_4={\displaystyle {\sum }_{i=1}^N{a}_i{P}_i}$ ; $\left(Q+Q_1+2\upsilon \right)=k^2\left(d_1+d_2\right)-2k\left(d_3+d_4{x}_{N+1}\right)+\left(x_{N+1}^2+T_{00}+2\upsilon \right)$ .

Thus,

$P\left(x_{N+1},k,\alpha ,\theta ,\varphi /x_1,x_2,\mathrm{...},x_n\right)\propto \exp \left(-\beta \left(\alpha -1\right)\right)C_1{\left[1+E_1{\left(k-C_2\right)}^2\right]}^{-d}$ (16)

where $C_1={\left\{x_{N+1}^2+T_{00}+2\upsilon -\left[{\left(d_3+d_4{x}_{N+1}\right)}^2/\left(d_1+d_2\right)\right]\right\}}^{}$ , $C_2=\left(d_3+d_4{x}_{N+1}\right)/\left(d_1+d_2\right)$ , $E_1=\left(d_1+d_2\right)/C_1$ .

Further, integrating out < $k$  from (16) we get

$P\left(x_{N+1},k,\alpha ,\theta ,\varphi /x_1,x_2,\mathrm{...},x_N\right)\propto \exp \left(-\beta \left(\alpha -1\right)\right)C_1^{-d}{E}_1^{-\left(1/2\right)}$                                   (17)

with $d=\left(\upsilon +1\right)/2$ which is the conditional predictive distribution of $x_{N+1}$  given $\alpha$ , $\theta$  and $\varphi$ . Further elimination of the parameters $\alpha$ , $\theta$  and $\varphi$ from (17) is not possible analytically. So the marginal posterior density of $x_{N+1}$  cannot be expressed in a closed form. Since the distribution in (17) is analytically not tractable, a complete Bayesian analysis is possible by numerical integration technique or simulation based approach, viz., MCMC technique.

Suppose one wants a point estimate (posterior mean) of $x_{N+1}$ , then one should compute the marginal posterior density of $x_{N+1}$  from (17) and use it to calculate the marginal posterior mean of $x_{N+1}$ . Thus four dimensional numerical integration is necessary in order to estimate $x_{N+1}$ . But it is a very difficult problem.

Practically, to perform four dimensional numerical integration is very difficult and therefore to reduce the dimensions of the numerical integration one may substitute the estimators, posterior means, $\widehat{\alpha }$ , $\widehat{\theta }$  and $\widehat{\varphi }$  respectively in the place of $\alpha$ , $\theta$  and $\varphi$ and then perform one dimensional numerical integration to find the conditional mean of $X_{N+1}$ . That is, one may eliminate the parameters as much as possible by analytical methods and then use the conditional estimates for the remaining parameters to compute the marginal posterior mean of the future observation.

A numerical example is considered for illustrating the one-step ahead predictive analysis of a future observation from the Canadian Lynx data. This data consists of the annual record of the numbers of Canadian Lynx trapped in the Mackenzie River district of North-west Canada for the period 1821 – 1934 (both years inclusive) giving a total of 114 observations. Brockwel and Davis⁵ (page 501) have transformed these data using the log transformation for the purpose of statistical analysis. These transformed data are used in our Bayesian predictive analysis.

Bayesian predictive distribution of the ${\left(r+1\right)}^{th}$ observation, using the r observation, is obtained. The mean of this distribution is taken to be the ${\left(r+1\right)}^{th}$  predicted value of the Lynx data. Since the direct evaluation of the mean of the one-step ahead predictive distribution involves four dimensional numerical integration, instead of the marginal predictive distribution of $X_{N+1}$ , the conditional predictive distribution of $X_{N+1}$ , given by (17) got by fixing the parameters $k$ , $\alpha$ , $\theta$  and $\varphi$  at their estimates, is used and the mean (posterior mean) is calculated using FORTRAN language. The posterior mean of the predictive distribution is computed numerically after fixing the parameters at their respective estimated values. The prediction is done for the cases r=11, 12,….,114 by taking first 10 observations as initial observations to estimate the parameters of the model and are given in the Table 1 which contains both the true values and the one-step ahead predicted values for the transformed data and the figure 1 represent graphically, the original data and one step-step ahead predicted values of the same. Figure 2 represent graphically, the original data for the last 14 observations and predicted values of the same through different methods, using FORTRAN program.