Diagnostic Analytics for an Autoregressive Model under the Skew-Normal Distribution

Abstract

Autoregressive models have played an important role in time series. In this paper, an autoregressive model based on the skew-normal distribution is considered. The estimation of its parameters is carried out by using the expectation–maximization algorithm, whereas the diagnostic analytics are conducted by means of the local influence method. Normal curvatures for the model under four perturbation schemes are established. Simulation studies are conducted to evaluate the performance of the proposed procedure. In addition, an empirical example involving weekly financial return data are analyzed using the procedure with the proposed diagnostic analytics, which has improved the model fit.

1. Introduction

For time series data analysis, autoregressive (AR) modeling is an essential technique and has been applied in many areas including biology, chemistry, earth sciences, economics, education, engineering, finance, health, medicine, and physics; see, for example, [1,2] for recent accounts of time series modeling and applications. Issues related to estimation and testing for AR models are extensive and well-established; see [3,4] for related issues to statistical diagnostics which are of equal importance. Local influence diagnostics is the study of how relevant minor perturbations impact the fit of the model and the results of statistical inference. This has become a useful statistical methodology after [5] introduced the idea of local influence to aid in the identification of potentially influential observations. Diagnostic analytics is used in a number of regression and time series models. Among others, Refs. [6,7,8,9,10,11,12,13,14,15] investigated the local influence of linear or nonlinear regression models under non-normal distributional assumptions. In a framework of time series data, Refs. [16,17,18] considered influence diagnostics for AR and vector AR models under normal or elliptical distributions.

A standard assumption for time series models is that their errors are mutually independent and follow normal or symmetric distributions, as studied in [16,17,18]. However, it is known that certain financial and other datasets feature errors with skewed distributions. In order to deal with such data, the skew-normal (SN) distribution and its scale-mixtures have provided an appealing alternative and can therefore be adopted. Their properties, extensions, and applications are becoming increasingly popular; see [19,20,21,22,23,24,25,26]. In addition, Refs. [27,28,29] studied SN linear, linear mixed, and nonlinear models. Ref. [8] analyzed diagnostics in the nonlinear model with scale mixtures of SN and AR errors. For particular financial applications, Ref. [30] looked at asset pricing issues with return data following SN models. To our knowledge, no study on influence diagnostics in AR time series models under SN distributions have been reported. Therefore, the objective of the present paper is to formulate an AR model under the SN distribution (SN AR model) and to derive diagnostic analytics with applications to financial data. We use the matrix differential calculus pioneered by [31] to establish the mathematical results used in our data analysis. We implement the maximum likelihood (ML) method with the expectation–maximization (EM) algorithm to estimate the SN AR model parameters, whereas the local influence method with four perturbation schemes is used for the diagnostic analytics. The EM algorithm has now become a popular iterative technique for the ML estimation method with incomplete data; see [32,33]. The paper proceeds as follows. In Section 2, the SN AR model is introduced and the ML estimations of the model parameters are derived. Section 3 presents the local influence method and establishes normal curvatures under four perturbation schemes. In Section 4, a simulation study and an empirical example involving an AR model are presented, while the effectiveness of the proposed diagnostics is illustrated and discussed. Our concluding remarks are addressed in Section 5. The derivations of the normal curvatures are presented in the Appendix A.

2. An SN AR Model and its Parameter Estimation

2.1. The SN AR Model

A random variable Y is said to follow an SN distribution with location parameter $\mu$, scale parameter ${\sigma }^2$ and skewness parameter $\lambda$, which is denoted by $Y\sim \mathrm{SN}(\mu ,{\sigma }^2,\lambda )$, if its probability density function is given by

$$f(y)=\frac{2}{\sigma }\varphi \left(\frac{y-\mu }{\sigma }\right)\Phi \left(\lambda \left(\frac{y-\mu }{\sigma }\right)\right),y\in \mathbb{R},\mu \in \mathbb{R},\sigma >0,\lambda \in \mathbb{R},$$

(1)

where $\varphi$ and $\Phi$ are the standard normal probability density and cumulative distribution functions, respectively. We see that, if $\lambda =0$, then the probability density function of Y defined in (1) reduces to the normal probability density function. If $Y\sim \mathrm{SN}(\mu ,{\sigma }^2,\lambda )$, then $\mathrm{E}(Y)=\mu +\sigma \delta \sqrt{2/\pi }$ and $\mathrm{Var}(Y)={\sigma }^2-(2/\pi ){\sigma }^2{\delta }^2$, where $\delta =\lambda /\sqrt{1+{\lambda }^2}$.

Let $Y_t$ be generated by a stationary AR(p) process given by

$$Y_t={\beta }_1{y}_{t-1}+\cdots +{\beta }_p{y}_{t-p}+u_t={\mathit{x}}_t^{\top }\mathit{\beta }+u_t,t=1,\dots ,T,$$

(2)

where $Y_t$ is a time series, with $Y_1,\dots ,Y_p$ being the p initial values for $y_{t-j}$, such that ${\mathit{x}}_t={(y_{t-1},\dots ,y_{t-p})}^{\top }$ is the $p\times 1$ vector of SN AR(p) covariates; ${\beta }_j$ is a regression coefficient, for $j=1,\dots ,p$, such that $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_p)}^{\top }$ is a $p\times 1$ regression coefficient vector, and $u_t$ is the scalar error term which follows an SN distribution, that is, $u_t\sim \mathrm{SN}(0,{\sigma }^2,\lambda )$, where ${\sigma }^2$ is the scale parameter and $\lambda$ is the skewness parameter. Thus, $\mathit{\theta }={(\mathit{\beta },{\sigma }^2,\lambda )}^{\top }$ is the $(p+2)\times 1$ vector of SN AR(p) model parameters. Note that we do not include an intercept in the SN model above. The expected value of $u_t$ is not zero, due to the assumption that the expected value of the underlying normal distribution is zero. Then, we choose the non-zero expected value of $u_t$ as an approximate to replace an intercept for the SN model.

2.2. ML Estimation

Finding the ML estimate of the parameter vector $\mathit{\theta }$ by direct maximization of the log-likelihood can potentially be a difficult task, so we implement the EM algorithm for this estimation. Let ${\mathit{y}}_c={({\mathit{y}}_o,{\mathit{y}}_m)}^{\top }$ denote the complete data, with ${\mathit{y}}_o$ being the observed data and ${\mathit{y}}_m$ the missing data. We use the notation ${\mathit{Y}}_c,{\mathit{Y}}_o,{\mathit{Y}}_m$ for the random vectors associated with ${\mathit{y}}_c$, ${\mathit{y}}_o$ and ${\mathit{y}}_m$, respectively. Starting from ${\mathit{\theta }}^{(0)}$ as the initial estimate, we can get ${\mathit{\theta }}^{(0)},{\mathit{\theta }}^{(1)},\dots$ by running the two steps of the EM algorithm iteratively as defined below.

In order to implement the EM algorithm, consider the stochastic representation given by

$$Y=\mu +\sigma \delta H+\sigma \sqrt{(1-{\delta }^2)}{H}_1,$$

(3)

where $H=|H_0|\sim \mathrm{HN}(0,1)$, with HN being used to indicate the half normal distribution. In addition, both $H_0$ and $H_1$ are independent random variables which follow a standard normal distribution. Note that, from (3), we have $Y|H=h\sim \mathrm{N}(\mu +h\lambda \sigma /\sqrt{1+{\lambda }^2},{\sigma }^2/(1+{\lambda }^2))$ and $H\sim \mathrm{HN}(0,1)$. Hence, by considering ${\mathit{y}}_o={(y_{p+1},\dots ,y_T)}^{\top }$, ${\mathit{y}}_m={(h_{p+1},\dots ,h_T)}^{\top }$ and ${\mathit{y}}_c={({\mathit{y}}_o,{\mathit{y}}_m)}^{\top }$ as the observed, missing and complete data sets, respectively, we have the complete-data log-likelihood function for $\mathit{\theta }={(\mathit{\beta },{\sigma }^2,\lambda )}^{\top }$ given by

$$\begin{array}{cc}\hfill {\ell }_c(\mathit{\theta },{\mathit{y}}_c)=& \sum _{t=p+1}^T\left(-\frac{1}{2}log({\sigma }^2)+\frac{1}{2}log(1+{\lambda }^2)-\frac{1+{\lambda }^2}{2{\sigma }^2}{\left(y_t-{\mathit{x}}_t^{\top }\mathit{\beta }-\frac{\lambda \sigma }{\sqrt{1+{\lambda }^2}}{h}_t\right)}^2\right).\hfill \end{array}$$

(4)

Therefore, for the E-step of the EM algorithm, given the current estimate ${\widehat{\mathit{\theta }}}^{(k)}$ and based on (4), we can calculate the function ${Q(\mathit{\theta })|}_{\mathit{\theta }={\widehat{\mathit{\theta }}}^{(k)}}=\mathrm{E}\left[{\ell }_c(\mathit{\theta },{\mathit{Y}}_c)\right|{\mathit{Y}}_o={\mathit{y}}_o]{|}_{\mathit{\theta }={\widehat{\mathit{\theta }}}^{(k)}}$ as

$$Q(\mathit{\theta })=-\frac{(T-p)}{2}log({\sigma }^2)+\frac{(T-p)}{2}log(1+{\lambda }^2)-\frac{(1+{\lambda }^2)}{2}{\sum _{t=p+1}^{\top }\left(\frac{y_t-{\mathit{x}}_t^{\top }\mathit{\beta }}{\sigma }-\frac{\lambda }{\sqrt{1+{\lambda }^2}}{\widehat{c}}_t\right)}^2-\frac{{\lambda }^2}{2}\sum _{t=p+1}^{\top }(\widehat{c_t^2}-{({\widehat{c}}_t)}^2),$$

(5)

where ${\widehat{c}}_t=\mathrm{E}(H_t|{\mathit{Y}}_o={\mathit{y}}_o){|}_{\mathit{\theta }={\widehat{\mathit{\theta }}}^{(k)}}={\tau }_1+{\tau }_2\varphi ({\tau }_1/{\tau }_2)/\Phi ({\tau }_1/{\tau }_2)$, $\widehat{c_t^2}=\mathrm{E}(H_t^2|{\mathit{Y}}_o={\mathit{y}}_o){|}_{\mathit{\theta }={\widehat{\mathit{\theta }}}^{(k)}}={\tau }_1^2+{\tau }_2^2+{\tau }_1{\tau }_2\varphi ({\tau }_1/{\tau }_2)/\Phi ({\tau }_1/{\tau }_2)$, and ${\tau }_1={({\widehat{\lambda }}^{(k)}/{\widehat{\sigma }}^{(k)})(1+{(\widehat{{\lambda }^{(k)}})}^2)}^{1/2}(y_t-{\mathit{x}}_t^{\top }{\widehat{\mathit{\beta }}}^{(k)})$, ${\tau }_2=1/{(1+{(\widehat{{\lambda }^{(k)}})}^2)}^{1/2}$.

Note that equations in (4) and (5) are obtained from the probability density function given in (1) and after calculating the expected value of (4). In addition, $\widehat{c_t^2}$ is different from ${({\widehat{c}}_t)}^2$. For the M-step, we update ${\widehat{\mathit{\theta }}}^{(k)}$ by the Newton–Raphson iteration $\dot{Q}({\widehat{\mathit{\theta }}}^{\left(k+1\right)})=\dot{Q}({\widehat{\mathit{\theta }}}^{(k)})+\ddot{Q}({\widehat{\mathit{\theta }}}^{(k)})({\widehat{\mathit{\theta }}}^{\left(k+1\right)}-{\widehat{\mathit{\theta }}}^{(k)})+o(|{\widehat{\mathit{\theta }}}^{\left(k+1\right)}-{\widehat{\mathit{\theta }}}^{(k)}|)$, where $\dot{Q}$ denotes the gradient vector and $\ddot{Q}$ denotes the Hessian matrix. As it is known, a suitable initial value ${\widehat{\mathit{\theta }}}^{(0)}$ is important and difficult to find in numerical computation. Thus, we can get ${\widehat{\mathit{\theta }}}^{(0)}=({\widehat{\mathit{\beta }}}^{(0)},{\widehat{{\sigma }^2}}^{(0)},{\widehat{\lambda }}^{(0)})$ considering ${\widehat{\mathit{\beta }}}^{(0)}$ as the ordinary least squares (OLS) estimate and then ${\widehat{{\sigma }^2}}^{(0)},{\widehat{\lambda }}^{(0)}$ can be computed as $\widehat{\mathit{\theta }}=(\widehat{\mathit{\beta }},\widehat{{\sigma }^2},\widehat{\lambda })$ until $|{\widehat{\mathit{\theta }}}^{\left(k+1\right)}-{\widehat{\mathit{\theta }}}^{(k)}|<{10}^{-5}$.

2.3. The Hessian Matrix

Next, we compute the Hessian matrix $\ddot{Q}\left(\mathit{\theta }\right)$ evaluated at the ML estimate $\widehat{\mathit{\theta }}$ using

$$\begin{array}{cc}\hfill {\ell }_c(\mathit{\theta },{\mathit{y}}_c)=& \sum _{t=p+1}^T\left(-\frac{1}{2}log({\sigma }^2)+\frac{1}{2}log(1+{\lambda }^2)-\frac{1+{\lambda }^2}{2{\sigma }^2}{\left(u_t-\frac{\lambda \sigma }{\sqrt{1+{\lambda }^2}}{h}_t\right)}^2\right).\hfill \end{array}$$

(6)

Then, based on (6), the Q function is given by

$$\begin{array}{c}{Q(\mathit{\theta })|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}=\mathrm{E}\left[{\ell }_c(\mathit{\theta },{\mathit{Y}}_c)\right|{\mathit{Y}}_o={\mathit{y}}_o]{|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}\hfill \\ =-\frac{(T-p)}{2}log({\sigma }^2)+\frac{(T-p)}{2}log(1+{\lambda }^2)\\ -\frac{(1+{\lambda }^2)}{2}{\sum _{t=p+1}^{\top }\left(\frac{u_t}{\sigma }-\frac{\lambda }{\sqrt{1+{\lambda }^2}}{\widehat{c}}_t\right)}^2-\frac{{\lambda }^2}{2}\sum _{t=p+1}^{\top }(\widehat{c_t^2}-{({\widehat{c}}_t)}^2).\hfill \end{array}$$

For the SN AR(p) model, we have the $(p+2)\times (p+2)$ Hessian matrix $\ddot{Q}(\widehat{\mathit{\theta }})$ given by

$$\ddot{Q}(\widehat{\mathit{\theta }})=\left(\begin{array}{ccc}\frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial \mathit{\beta }\partial {\mathit{\beta }}^{\top }}& \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial \mathit{\beta }\partial {\sigma }^2}& \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial \mathit{\beta }\partial \lambda }\\ \\ \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial {\sigma }^2\partial \mathit{\beta }}& \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial {\left({\sigma }^2\right)}^2}& \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial {\sigma }^2\partial \lambda }\\ \\ \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial \lambda \partial \mathit{\beta }}& \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial \lambda \partial {\sigma }^2}& \frac{{\partial }^{2}Q\left(\mathit{\theta }\right)}{\partial {\lambda }^2}\end{array}\right){|}_{\mathit{\theta }=\widehat{\mathit{\theta }}},$$

(7)

where $\widehat{\mathit{\theta }}=(\widehat{\mathit{\beta }},{\widehat{\sigma }}^2,\widehat{\lambda })$ is the ML estimate of $\mathit{\theta }$. The expression for ${\partial }^{2}Q\left(\mathit{\theta }\right)/\partial \mathit{\beta }\partial {\mathit{\beta }}^{\top }$ and the other submatrices are presented in Appendix A.1.

3. Influence Analysis for the SN AR Model

3.1. Local Influence

Let $\ell (\mathit{\theta })$ denote the log-likelihood function for the model given in (2), which is the postulated model, where $\mathit{\theta }$ is a $(p+2)\times 1$ vector of unknown parameters with its ML estimate $\widehat{\mathit{\theta }}$. Let $\mathit{\omega }={({\omega }_1,\dots ,{\omega }_q)}^{\top }$ denote a $q\times 1$ vector of perturbations confined to some open subset of ${\mathbb{R}}^q$ and let ${\mathit{\omega }}_0$ denote a $q\times 1$ non-perturbation vector, with q as a suitable number of dimensions and ${\mathit{\omega }}_0=(0,\dots ,0)$, or ${\mathit{\omega }}_0=(1,\dots ,1)$, or a third choice, depending on the context. Then, $\ell (\mathit{\theta })$ and $\ell (\mathit{\theta }|\mathit{\omega })$ represent the log-likelihood functions of the postulated and perturbed models, respectively. Note that $\ell (\mathit{\theta })=\ell (\mathit{\theta }|{\mathit{\omega }}_0)$.

We suppose that $\ell (\mathit{\theta }|\mathit{\omega })$ is twice continuously differentiable in a neighborhood of $(\widehat{\mathit{\theta }},{\mathit{\omega }}_0)$. We are interested in the comparison of $\widehat{\mathit{\theta }}$ and ${\widehat{\mathit{\theta }}}_{\mathit{\omega }}$ using the local influence method idea, which is to investigate the degree of inference affected by those minor changes in the corresponding perturbations. Ref. [5] used likelihood displacement (LD) to assess the influence of the perturbation $\mathit{\omega }$ defined as $\mathrm{LD}(\mathit{\omega })=2(\ell (\widehat{\mathit{\theta }})-\ell ({\widehat{\mathit{\theta }}}_{\mathit{\omega }}))$. Here, large values of $\mathrm{LD}(\mathit{\omega })$ indicate that $\widehat{\mathit{\theta }}$ and ${\widehat{\mathit{\theta }}}_{\mathit{\omega }}$ differ considerably in relation to the contours of the non-perturbed log-likelihood function $\ell (\mathit{\theta })$. The idea is based on studying the local behavior of $\mathrm{LD}(\mathit{\omega })$ and the normal curvature $C_l\left(\mathit{\theta }\right)$ in a unit-length direction vector l, where $\left|\right|l\left|\right|=1$. According to [5], the normal curvature used to examine the local influence of the perturbation vector at $\mathit{\omega }={\mathit{\omega }}_0$ is $C_l\left(\mathit{\theta }\right)=2|l^{\top }\ddot{F}l|=2|l^{\top }({\Delta }^{\top }{\ddot{\ell }}^{-1}\Delta )l|$, with

$$\ddot{F}=\frac{{\partial }^2\ell (\mathit{\theta }|\mathit{\omega })}{\partial \mathit{\omega }\partial {\mathit{\omega }}^{\top }},\Delta =\frac{{\partial }^2\ell (\mathit{\theta }|\mathit{\omega })}{\partial \mathit{\theta }\partial {\mathit{\omega }}^{\top }},\ddot{\ell }=\frac{{\partial }^2\ell (\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }},$$

where l is a $q\times 1$ vector of unit length, $-\ddot{\ell }$ is the $(p+2)\times (p+2)$ observed information matrix for the postulated model, $\Delta$ is the $(p+2)\times q$ perturbation matrix for the perturbed model, and $-\ddot{\ell }$ and $\Delta$ are evaluated at $\mathit{\theta }=\widehat{\mathit{\theta }}$ and $\mathit{\omega }={\mathit{\omega }}_0$. The suggestion is to make the local influence diagnostic analytics by finding the maximum curvature $C_{max}={max}_{\left|\right|l\left|\right|=1}{C}_l$, where $C_{max}$ corresponds to the largest absolute eigenvalue ${\lambda }_{max}$ and its associated eigenvector $l_{max}$ of the matrix $\ddot{F}={\Delta }^{\top }{\ddot{\ell }}^{-1}\Delta$. If the absolute value of the ith element of $l_{max}$ is the largest, then the ith observation in the data may be the most influential. To examine the magnitude of influence, it is useful to have a benchmark value for $C_{max}$ and for the elements of $l_{max}$, see [10,18,34].

3.2. Local Influence for the SN AR Model

Next, we conduct a local influence diagnostic analytics for the SN AR(p) model. Due to the complexity of the SN distribution, we obtain the ML estimates based on the EM algorithm. As suggested by [29,34], the Q function and a Q function obtained similarly as LD can be used to replacethe log-likelihood function and likelihood displacement in the method of [5], in order to assess the influence of the perturbation. Thus, the normal curvature should be changed to be

$$C_l\left(\mathit{\theta }\right)=2|l^{\top }\ddot{F}l|=2|l^{\top }({\Delta }^{\top }{\ddot{Q}}^{-1}\Delta )l|,$$

(8)

with

$$\ddot{F}=\frac{{\partial }^{2}Q(\mathit{\theta }|\mathit{\omega })}{\partial \mathit{\omega }\partial {\mathit{\omega }}^{\top }},\Delta =\frac{{\partial }^{2}Q(\mathit{\theta }|\mathit{\omega })}{\partial \mathit{\theta }\partial {\mathit{\omega }}^{\top }},\ddot{Q}=\frac{{\partial }^{2}Q(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }},$$

where l is a $q\times 1$ vector of unit length, and $\ddot{F}$, $\ddot{Q}$ and $\Delta$ are $q\times q$, $(p+2)\times (p+2)$ and $(p+2)\times q$ matrices, respectively. In addition, $\ddot{Q}$ and $\Delta$ need to be evaluated at $\mathit{\theta }=\widehat{\mathit{\theta }}$ and $\mathit{\omega }={\mathit{\omega }}_0$.

The method examines the total local influence $C_t=C_{l_t}(\mathit{\theta })$, where $l_t$ is a $q\times 1$ unit-length vector with one at the tth position and zeros elsewhere. We denote $S=-{\Delta }^{\top }{\ddot{Q}}^{-1}\Delta$. Since $C_l\left(\mathit{\theta }\right)$ is not invariant under a uniform change of scale, Ref. [34] proposed the conformal normal curvature $B_l(\mathit{\theta })=C_l(\mathit{\theta })/(2\mathrm{tr}(S)$). An interesting property of the conformal normal curvature is that, for any unit-length direction l, we have $0\le B_l\left(\mathit{\theta }\right)\le 1$, which allows comparison of curvatures among different models.

In order to determine if the tth observation is influential, Ref. [34] proposed to classify the tth observation as possibly influential if $M{(0)}_t=B_{l_t}$ is greater than the benchmark $1/q+c^{\ast }\mathrm{SM}(0)$, where $\mathrm{SM}(0)$ is the sample standard error of $M{(0)}_k$, for $k=1,\dots ,q$, and $c^{\ast }$ is a certain constant. Depending on the specific application, $c^{\ast }$ may be taken to be a suitably selected positive value.

The forms given in (8) are used to obtain our normal curvature results under the four perturbation schemes, namely the case-weights, data, variance parameter, and skewness parameter schemes. The matrices $\ddot{Q}$ and $\Delta$ need to be established for each scheme.

We employ the matrix differential calculus proposed by [31] to establish these algebraic results in the following sections. We present their derivations in their respective Appendix A.

3.3. Perturbation Matrices

3.3.1. Perturbation of Case-Weights

Assume that a minor perturbation is made on the SN AR(p) model, with $y_t={\mathit{x}}_t^{\top }\mathit{\beta }+u_t$ being replaced by ${\omega }_t{y}_t={\omega }_t{\mathit{x}}_t^{\top }\mathit{\beta }+{\omega }_t{u}_t$, where ${\omega }_t$ is the weight. Let $\mathit{\omega }={({\omega }_{p+1},\dots ,{\omega }_T)}^{\top }$ denote the $(T-p)\times 1$ perturbation vector and ${\mathit{\omega }}_0={(1,\dots ,1)}^{\top }$ denote the $(T-p)\times 1$ non-perturbation vector. Then, the complete-data log-likelihood function of the perturbed model is given by

$${\ell }_c(\mathit{\theta },\mathit{\omega },{\mathit{y}}_c)=\sum _{t=p+1}^T\left(-\frac{1}{2}log({\sigma }^2)+\frac{1}{2}log(1+{\lambda }^2)-\frac{1+{\lambda }^2}{2{\sigma }^2}{\left({\omega }_t{u}_t-\frac{\lambda \sigma }{\sqrt{1+{\lambda }^2}}{h}_t\right)}^2\right).$$

(9)

Thus, the perturbed Q function obtained from (9) is expressed as

$$\begin{array}{c}{Q(\mathit{\theta },\mathit{\omega })|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}=\mathrm{E}\left[{\ell }_c(\mathit{\theta },\mathit{\omega },{\mathit{Y}}_c)\right|{\mathit{Y}}_o={\mathit{y}}_o]{|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}\hfill \\ =-\frac{(T-p)}{2}log({\sigma }^2)+\frac{(T-p)}{2}log(1+{\lambda }^2)\\ \hfill -\frac{1+{\lambda }^2}{2}{\sum _{t=p+1}^{\top }\left(\frac{{\omega }_t{u}_t}{\sigma }-\frac{\lambda }{\sqrt{1+{\lambda }^2}}{\widehat{c}}_t\right)}^2-\frac{{\lambda }^2}{2}\sum _{t=p+1}^{\top }(\widehat{c_t^2}-{({\widehat{c}}_t)}^2).\end{array}$$

(10)

For the SN AR$(p)$ model in the perturbation of case-weights, the $(p+2)\times (T-p)$ perturbation matrix $\Delta$ must be evaluated at $\mathit{\theta }=\widehat{\mathit{\theta }}$ and $\mathit{\omega }={\mathit{\omega }}_0$ after taking derivatives, obtaining

$$\Delta =\frac{{\partial }^{2}Q(\mathit{\theta },\mathit{\omega })}{\partial \mathit{\theta }\partial \mathit{\omega }}{|}_{\mathit{\theta }=\widehat{\mathit{\theta }},\mathit{\omega }={\mathit{\omega }}_0}=\left(\begin{array}{c}\frac{1+{\widehat{\lambda }}^2}{{\widehat{\sigma }}^2}{\widehat{u}}_t{{\mathit{x}}_t}^{\top }\\ \frac{1+{\widehat{\lambda }}^2}{2{\widehat{\sigma }}^3}\left(\frac{2{\widehat{u}}_t^2}{\widehat{\sigma }}-\frac{\widehat{\lambda }{\widehat{c}}_t}{\sqrt{1+{\widehat{\lambda }}^2}}{\widehat{u}}_t\right)\\ -\frac{2\widehat{\lambda }}{{\widehat{\sigma }}^2}{\widehat{u}}_t^2+\frac{2\widehat{\lambda }+1}{\sqrt{1+{\widehat{\lambda }}^2}}\frac{{\widehat{c}}_t}{\widehat{\sigma }}{\widehat{u}}_t\end{array}\right),$$

(11)

where

$$\begin{array}{cc}\hfill & {\widehat{u}}_t=y_t-{\mathit{x}}_t^{\top }\widehat{\mathit{\beta }},\hfill \\ \hfill & {\widehat{c}}_t=\mathrm{E}(H_t|Y_0,\widehat{\mathit{\theta }})={\tau }_1+\frac{\varphi (\frac{{\tau }_1}{{\tau }_2})}{\Phi (\frac{{\tau }_1}{{\tau }_2})}{\tau }_2,\hfill \\ \hfill & \widehat{c_t^2}=\mathrm{E}({H_t}^2|Y_0,\widehat{\mathit{\theta }})={\tau }_1^2+{\tau }_2^2+\frac{\varphi (\frac{{\tau }_1}{{\tau }_2})}{\Phi (\frac{{\tau }_1}{{\tau }_2})}{\tau }_1{\tau }_2,\hfill \\ \hfill & {\tau }_1=\frac{\widehat{\lambda }}{\widehat{\sigma }\sqrt{1+{\widehat{\lambda }}^2}}{\widehat{u}}_t,{\tau }_2=\frac{1}{\sqrt{1+{\widehat{\lambda }}^2}}.\hfill \end{array}$$

3.3.2. Perturbation of Data

Assume that a perturbation replaces $y_t$ by ${\omega }_t+y_t$. Let $\mathit{\omega }={({\omega }_{p+1},\dots ,{\omega }_T)}^{\top }$ denote the $(T-p)\times 1$ perturbation vector and ${\mathit{\omega }}_0={(0,\dots ,0)}^{\top }$ the $(T-p)\times 1$ non-perturbation vector. The perturbed AR$(p)$ model can be written as $y_t+{\omega }_t={\beta }_1\left(y_{t-1}+{\omega }_{t-1}\right)+\cdots +{\beta }_p\left(y_{t-p}+{\omega }_{t-p}\right)+u_t$, where $u_t=y_t-{\mathit{x}}_t^{\top }\mathit{\beta }+\mu ({\omega }_t)$ and $\mu ({\omega }_t)={\omega }_t-{\beta }_1{\omega }_{t-1}-\cdots -{\beta }_p{\omega }_{t-p}$. Then, the complete-data log-likelihood function of the perturbed model is given by

$${\ell }_c(\mathit{\theta },\mathit{\omega },{\mathit{y}}_c)=\sum _{t=p+1}^T\left(-\frac{1}{2}log({\sigma }^2)+\frac{1}{2}log(1+{\lambda }^2)-\frac{1+{\lambda }^2}{2{\sigma }^2}{\left(u_t+\mu ({\omega }_t)-\frac{\lambda \sigma }{\sqrt{1+{\lambda }^2}}{h}_t\right)}^2\right).$$

Thus, the perturbed Q function is expressed as

$$\begin{array}{ccc}\hfill {Q(\mathit{\theta },\mathit{\omega })|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}& =\hfill & -\frac{(T-p)}{2}log({\sigma }^2)+\frac{(T-p)}{2}log(1+{\lambda }^2)\hfill \\ & & -\frac{1+{\lambda }^2}{2}{\sum _{t=p+1}^{\top }\left(\frac{u_t+\mu ({\omega }_t)}{\sigma }-\frac{\lambda }{\sqrt{1+{\lambda }^2}}{\widehat{c}}_t\right)}^2-\frac{{\lambda }^2}{2}\sum _{t=p+1}^{\top }(\widehat{c_t^2}-{({\widehat{c}}_t)}^2).\hfill \end{array}$$

(12)

For the SN AR$(p)$ model in the perturbation of data, we have

$$\Delta =\frac{{\partial }^{2}Q(\mathit{\theta },\mathit{\omega })}{\partial \mathit{\theta }\partial \mathit{\omega }}{|}_{\mathit{\theta }=\widehat{\mathit{\theta }},\mathit{\omega }={\mathit{\omega }}_0}=\left(\begin{array}{c}\frac{1+{\widehat{\lambda }}^2}{{\widehat{\sigma }}^2}{{\mathit{x}}_t}^{\top }\\ \frac{1+{\widehat{\lambda }}^2}{2{\widehat{\sigma }}^3}\left(\frac{2{\widehat{u}}_t^2}{\widehat{\sigma }}-\frac{\widehat{\lambda }}{\sqrt{1+{\widehat{\lambda }}^2}}{\widehat{c}}_t\right)\\ -\frac{2\widehat{\lambda }}{{\widehat{\sigma }}^2}{\widehat{u}}_t^2+\frac{2\widehat{\lambda }+1}{\sqrt{1+{\widehat{\lambda }}^2}}\frac{{\widehat{c}}_t}{\widehat{\sigma }}\end{array}\right),$$

(13)

where ${\widehat{u}}_t=y_t-{\mathit{x}}_t^{\top }\widehat{\mathit{\beta }}$, and ${\widehat{h}}_t$ is as in the data perturbation.

3.3.3. Perturbation of Scale

Assume that the variance parameter ${\sigma }^2$ in the model is replaced by ${{\omega }_t}^{-1}{\sigma }^2$, meaning that $u_t\sim \mathrm{SN}(0,{{\omega }_t}^{-1}{\sigma }^2,\lambda )$. Let $\mathit{\omega }={({\omega }_{p+1},\dots ,{\omega }_T)}^{\top }$ denote the $(T-p)\times 1$ perturbation vector and ${\mathit{\omega }}_0={(1,\dots ,1)}^{\top }$ denote the $(T-p)\times 1$ non-perturbation vector. Then, the complete-data log-likelihood function of the perturbed model is given by

$${\ell }_c(\mathit{\theta },\mathit{\omega },{\mathit{y}}_c)=\sum _{t=p+1}^T\left(-\frac{1}{2}log({\sigma }^2)+\frac{1}{2}log({\omega }_t)+\frac{1}{2}log(1+{\lambda }^2)-\frac{{\omega }_t(1+{\lambda }^2)}{2{\sigma }^2}{\left(u_t-\frac{\lambda \sigma }{\sqrt{{\omega }_t(1+{\lambda }^2)}}{h}_t\right)}^2\right).$$

Thus, the perturbed Q function is expressed as

$$\begin{array}{ccc}\hfill {Q(\mathit{\theta },\mathit{\omega })|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}& =\hfill & -\frac{(T-p)}{2}log({\sigma }^2)+\sum _{t=p+1}^{\top }\frac{1}{2}log({\omega }_t)+\frac{(T-p)}{2}log(1+{\lambda }^2)\hfill \\ & & -\frac{1+{\lambda }^2}{2}{\sum _{t=p+1}^{\top }\left(\frac{\sqrt{{\omega }_t}{u}_t}{\sigma }-\frac{\lambda }{\sqrt{1+{\lambda }^2}}{\widehat{c}}_t\right)}^2-\frac{{\lambda }^2}{2}\sum _{t=p+1}^{\top }(\widehat{c_t^2}-{({\widehat{c}}_t)}^2).\hfill \end{array}$$

(14)

For the SN AR$(p)$ model with the perturbation of variance parameter, we have

$$\Delta =\frac{{\partial }^{2}Q(\mathit{\theta },\mathit{\omega })}{\partial \mathit{\theta }\partial \mathit{\omega }}{|}_{\mathit{\theta }=\widehat{\mathit{\theta }},\mathit{\omega }={\mathit{\omega }}_0}=\left(\begin{array}{c}\frac{1+{\widehat{\lambda }}^2}{\widehat{\sigma }}\left(\frac{\widehat{u}{}_t}{\widehat{\sigma }}-\frac{\widehat{\lambda }}{2\sqrt{1+{\widehat{\lambda }}^2}}{\widehat{c}}_t\right){{\mathit{x}}_t}^{\top }\\ \frac{\left(1+{\widehat{\lambda }}^2\right){\widehat{u}}_t^2}{2{\widehat{\sigma }}^4}-\frac{\widehat{\lambda }\sqrt{1+{\widehat{\lambda }}^2}}{4{\widehat{\sigma }}^3}{\widehat{c}}_t{\widehat{u}}_t\\ -\frac{\widehat{\lambda }}{{\widehat{\sigma }}^2}{\widehat{u}}_t^2+\frac{2{\widehat{\lambda }}^2+1}{\sqrt{1+{\widehat{\lambda }}^2}}\frac{{\widehat{c}}_t}{2\widehat{\sigma }}{\widehat{u}}_t\end{array}\right),$$

(15)

where ${\widehat{u}}_t=y_t-{\mathit{x}}_t^{\top }\widehat{\mathit{\beta }}$, and ${\widehat{h}}_t$ is as in data perturbation.

3.3.4. Perturbation of Skewness

Considering the particular skewed feature of the distribution, we may investigate the effect on the model fit by making a minor change of the skewness parameter $\lambda$. In our perturbed model, we propose to replace $\lambda$ by ${\sqrt{\omega }}_i\lambda$. Let $\mathit{\omega }={({\omega }_{p+1},\dots ,{\omega }_T)}^{\top }$ denote the $(T-p)\times 1$ perturbation vector and ${\omega }_0={(1,\dots ,1)}^{\top }$ denote the $(T-p)\times 1$ non-perturbation vector. Then, the complete-data log-likelihood function of the perturbed model is given by

$${\ell }_c(\mathit{\theta },\mathit{\omega },{\mathit{y}}_c)=\sum _{t=p+1}^T\left(-\frac{1}{2}log({\sigma }^2)+\frac{1}{2}log(1+{{\omega }_t\lambda }^2)-\frac{1+{{\omega }_t\lambda }^2}{2{\sigma }^2}{\left(u_t-\frac{\sqrt{{\omega }_t}\lambda \sigma }{\sqrt{1+{{\omega }_t\lambda }^2}}{h}_t\right)}^2\right).$$

Thus, the perturbed Q function is expressed as

$$\begin{array}{ccc}\hfill {Q(\mathit{\theta },\mathit{\omega })|}_{\mathit{\theta }=\widehat{\mathit{\theta }}}& =\hfill & -\frac{(T-p)}{2}log({\sigma }^2)+\frac{(T-p)}{2}log(1+{\omega }_t{\lambda }^2)\hfill \\ & & -{\sum _{t=p+1}^{\top }\frac{1+{\omega }_t{\lambda }^2}{2}\left(\frac{u_t}{\sigma }-\frac{\sqrt{{\omega }_t}\lambda }{\sqrt{1+{\omega }_t{\lambda }^2}}{\widehat{c}}_t\right)}^2-\sum _{t=p+1}^{\top }\frac{{\omega }_t{\lambda }^2}{2}(\widehat{c_t^2}-{({\widehat{c}}_t)}^2).\hfill \end{array}$$

For the SN AR$(p)$ model with the perturbation of the skewness parameter, we have

$$\Delta =\frac{{\partial }^{2}Q(\mathit{\theta },\mathit{\omega })}{\partial \mathit{\theta }\partial \mathit{\omega }}{|}_{\mathit{\theta }=\widehat{\mathit{\theta }},\mathit{\omega }={\mathit{\omega }}_0}=\left(\begin{array}{c}\left(\frac{{\widehat{u}}_t}{{\widehat{\sigma }}^2}{\widehat{\lambda }}^2-\frac{\widehat{\lambda }\left(2{\widehat{\lambda }}^2+1\right)}{2\sqrt{1+{\widehat{\lambda }}^2}}\frac{{\widehat{c}}_t}{\widehat{\sigma }}\right){{\mathit{x}}_t}^{\top }\\ \frac{{\widehat{\lambda }}^2{\widehat{u}}_t^2}{2{\widehat{\sigma }}^4}-\frac{\widehat{\lambda }}{2{\widehat{\sigma }}^3}\frac{1+2{\widehat{\lambda }}^2}{2\sqrt{1+{\widehat{\lambda }}^2}}{\widehat{c}}_t{\widehat{u}}_t\\ \frac{\widehat{\lambda }}{{\left({\widehat{\lambda }}^2+1\right)}^2}-\frac{\widehat{\lambda }}{{\widehat{\sigma }}^2}{\widehat{u}}_t^2-\widehat{\lambda }{{\widehat{c}}_t}^2+\frac{\left(6{\widehat{\lambda }}^2+1\right)\left({\widehat{\lambda }}^2+1\right)-{\widehat{\lambda }}^2\left(2{\widehat{\lambda }}^2+1\right)}{2{\left({\widehat{\lambda }}^2+1\right)}^{3/2}}\frac{{\widehat{c}}_t}{\widehat{\sigma }}{\widehat{u}}_t\end{array}\right),$$

(16)

where ${\widehat{u}}_t=y_t-{\mathit{x}}_t^{\top }\widehat{\mathit{\beta }}$ and ${\widehat{h}}_t$, $\widehat{c_t^2}$ are as in the data perturbation.

4. Numerical Results

4.1. Simulation Study I: Effectiveness of the Diagnostics

For our simulation, we consider an AR$(1)$ model as specified in Section 2 given by

$$y_t=\beta y_{t-1}+u_t,u_t\sim \mathrm{SN}(0,{\sigma }^2,\lambda ),$$

where we choose $\beta =0.12$, ${\sigma }^2=0.003$ and $\lambda =0.1$. We generate $T=400$ observations. Now, we compare the performance of the ML estimates in the presence of five perturbed or shifted observations among three different scenarios with $\lambda =0.1,0.2,0.3$.

We perturb the value $y_t$ by ${y_t}^{\ast }=y_t+\beta y_{t-1}d$, where $t=200,201,202,203,204$ and $d=5,10,\dots ,50$ in order to guarantee the presence of outliers. We calculate the ML estimate of $\beta$ under the perturbed and non-perturbed data, by fitting both data sets under the AR(1) model with $\lambda =0.1,0.2,0.3$, respectively. Then, we compute the relative changes of the estimates as $|(E_{st(i)}-E_{st})/E_{st}|$, where $E_{st(i)}$ is the estimate of $\beta$ under the perturbed data and $E_{st}$ is the estimate of $\beta$ under the non-perturbed data. The above computation in MATLAB (version 9.3, MathWorks, US) runs five iterations to converge in less than 60 s on a PC. Our experience indicates that the computation runs up to 10 iterations to converge in less than 120 s for even $T=1000$. Figure 1 shows the effectiveness of the influence diagnostics.

4.2. Simulation Study II: Comparing SN and Normal Distributions

We perform a numerical simulation to examine the effectiveness of our method under SN distributions. The results under the SN and normal distributions are compared as follows:

• Step 1. We use the simulated data $(\lambda =0.1)$ with $y_t$ perturbed by ${y_t}^{\ast }=y_t+\beta y_{t-1}d$, where $d=5$ and $t=200,201,202,203,204$. We fit an AR(1) model under normality to the data, and obtain the fitted AR(1) model given by $Y_t=0.1549y_{t-1}+u_t,$ where $u_t\sim \mathrm{N}(0,0.0151)$.
• Step 2. We conduct a local influence diagnostic analytics under the normal distribution using the diagnostic results given by [18].
• Step 3. We compare the local influence results for the normal distribution in Step 2. In Figure 2, twenty-four influence observations are detected under the SN distribution. These results are summarized in

  Table 1. Local influence results for normal and SN models.

  ID	Index under the Normal Model	Index under the SN Model	Observed Value
  1	33	33	−0.334
  2	34	34	0.030
  3	-	62	−0.271
  4	77	77	−0.260
  5	111	111	−0.297
  6	112	112	−0.225
  7	140	140	0.319
  8	141	141	0.089
  9	-	153	−0.175
  10	-	201	−0.269
  11	202	202	−0.326
  12	203	203	−0.090
  13	205	205	0.328
  14	206	206	−0.124
  15	214	214	−0.406
  16	215	215	−0.159
  17	293	293	0.296
  18	294	294	−0.012
  19	-	301	−0.278
  20	306	306	0.043
  21	345	345	−0.293
  22	346	346	−0.092
  23	362	362	−0.307
  24	363	363	−0.048

  .

As presented, twenty-four influential observations are detected by the local influence diagnostic analytics under the SN distribution, which is more than the twenty influential values detected by the local influence analytics under the normal distribution. These results are informative and as expected. We see that the influential values (that is, cases #62, #153, #201 and #301) are less than zero.

As we understand, when $\lambda$ for the SN distribution is larger than zero, it is easier to find possible influential values less than zero due to the difference in patterns between the SN distribution and the normal distribution. This indicates that the diagnostic results under the SN distribution established in Section 3 work well.

4.3. Real-World Data Analysis

We conduct an empirical example of financial data based on our results presented in Section 2 and Section 3. We choose Chevron shares (hereinafter referred to as CVX weekly financial return data), which were collected from 12 January 2007 to 1 August 2014 to construct the AR models. Figure 3 shows the log transformation of this data set (a total of 395 observations).

We first determine the order of the AR model, using the following steps:

• Step 1. We assume the data are subject to an AR$(p)$ model, for $p=1,2,\dots$, given by

  $$\begin{array}{c}{Y}_t={\beta }_1{y}_{t-1}+u_t\hfill \\ ⋮\hfill \\ Y_t={\beta }_1{y}_{t-1}+{\beta }_2{y}_{t-2}+\dots +{\beta }_p{y}_{t-p}+u_t\hfill \\ ⋮\hfill \end{array}$$

  (17)
• Step 2. For the ith equation given in (17), we let ${\widehat{\beta }}_j^{(i)}$ be the OLS estimate of ${\beta }_j$, where the superscript $(i)$ denotes the estimates for the AR(i) model. Then, the residual is defined as ${\widehat{u}}_t^{(i)}=y_t-{\widehat{\beta }}_1^{(i)}{y}_{t-1}-{\widehat{\beta }}_2^{(i)}{y}_{t-2}-\cdots -{\widehat{\beta }}_i^{(i)}{y}_{t-i}.$ The estimated variance for the AR(i) model is expressed as ${\widehat{\sigma }}_i^2=(1/T-2i-1){\sum }_{t=i+1}^{\top }{({\widehat{u}}_t^{(i)})}^2$.
• Step 3. We use the ith and $(i-1)$th equations given in (17) to test ${\mathrm{H}}_0:{\beta }_i=0$ versus ${\mathrm{H}}_1:{\beta }_i\ne 0$, that is, we test the AR(i) model versus the AR($i-1$) model. The test statistic is defined as

  $$M(i)=-(T-i-2.5)log\left(\frac{{{\widehat{\sigma }}_i}^2}{{{\widehat{\sigma }}_{i-1}}^2}\right).$$

  (18)

For our model, $M(i)$ follows an asymptotically chi-square distribution with one degree of freedom, that is, $M(i)\sim {\chi }^2(1)$.

We calculate $M(i)$ by (18), for $i=1,\dots ,7$, and present the result in

Table 2. Test statistic $M(i)$, for $i=1,\dots ,7$ using CVX data.

Order	1	2	3	4	5	6	7
$M(i)$	−0.3826	−1.001	16.9551	2.9495	−1.6675	−1.9619	8.9958

. As the 99th percentile of the chi-square distribution with one degree of freedom is 6.635 (${\chi }_{0.99}^2(1)=6.635$ from Table 2), we select the order of the AR(p) model to be $p=3$.

From the EM algorithm detailed in Section 2, we can find the estimate of parameter $\widehat{\mathit{\theta }}$ given by

$$(\widehat{{\beta }_1},\widehat{{\beta }_2},\widehat{{\beta }_3},{\widehat{\sigma }}^2,\widehat{\lambda })=(0.0153,-0.0333,-0.2198,0.0013,0.0661).$$

Since the absolute values of ${\beta }_1$, ${\beta }_2$ and ${\beta }_3$ all are less than 1, we can accept that the CVX time series is stationary for the AR(3) model. Thus, we obtain a predictive model in the fitted SN AR(3) structure given by

$${\widehat{y}}_t=0.0153y_{t-1}-0.0333y_{t-2}-0.2198y_{t-3},$$

with $\widehat{\mu }=0$, ${\widehat{\sigma }}^2=0.0013$, and $\widehat{\lambda }=0.0661$.

By applying the method in Section 3, we can conduct influence analysis for the SN AR(3) model. After calculating the observed Hessian matrix and the $\Delta$ matrices for the four perturbation schemes of case-weights, data, variance parameter, and skewness parameter, we obtain the diagnostics matrices $S_1$, $S_2$, $S_3$, and $S_4$, respectively. In this case, we select the benchmark as $1/392+3\mathrm{SM}(0)$, with the values of 0.1144, 0.0178, 0.0347, and 0.0791 established in the simulation study for the four perturbation schemes, respectively.

In Figure 4, the straight line represents the benchmark value which determines whether an observation is potentially influential. The observation is justified to be potentially influential when its diagnostic value exceeds the benchmark.

Firstly, we only find case #92 to be potentially influential. The other potential influential observations may be masked by case #92. Similar to a step-wise diagnostic procedure proposed by [6], a second step of identification of influential observations is conducted. In the second rstep, we replace the value of case #92 by the average of cases #91 and #93 to get a new time series.

We re-fit an AR model in the same manner as done previously in the first step. For the new time series, the parameters in the SN AR model are again estimated by applying the EM-algorithm with the order also selected to be three. Thus, we present the new SN AR(3) model as

$${\widehat{y}}_t=0.0092y_{t-1}-0.0374y_{t-2}-0.1328y_{t-3},$$

with $\widehat{\mu }=0$, ${\widehat{\sigma }}^2=0.0011$, $\widehat{\lambda }=0.1015$. Since the absolute values of ${\widehat{\beta }}_1{\widehat{\beta }}_2$, and ${\widehat{\beta }}_3$ are all less than one, we can accept the CVX time series to be stationary for the AR(3) model. For the new AR(3) model, we conduct influence analysis. We select the benchmark as $1/392+3\mathrm{SM}(0)$ again, with the four values of 0.0377, 0.0100, 0.0162, and 0.0203 for the corresponding perturbation schemes, respectively. Thus, twenty-two influence observations are detected in Figure 5, summarized in

Table 3. Summary of the curvature-based diagnostic analytics.

ID	Time	CVXs Return	Case-Weights	Data	Variance	Skewness
54	2008-01-18	−8.286%				*
79	2008-07-11	−6.687%				*
87	2008-09-05	−7.329%				*
91	2008-10-03	−9.109%				*
92	2008-10-10	−31.674%	*	*	*	*
95	2008-10-31	15.4665%	*		*
96	2008-11-07	−1.54%		*
97	2008-11-14	−1.067%		*
98	2008-11-21	−3.06%		*
99	2008-11-28	11.4104%	*
101	2008-12-12	5.9723%		*
102	2008-12-19	−10.888%	*	*	*	*
103	2008-12-26	−0.708%		*
104	2009-01-02	8.4069%			*
105	2009-01-09	−4.956%		*		*
110	2009-02-13	−7.152%				*
113	2009-03-06	−4.102%		*
181	2010-06-25	−7.505%				*
246	2011-09-23	−10.154%	*			*
254	2011-11-18	−8.955%				*
256	2011-12-02	9.6993%			*
257	2011-12-09	2.4863%		*

. It tallies the observations which are identified to be potentially influential for CVX, and an "*" in Table 3 indicates that the observation has been identified via the assigned perturbation scheme. It is worth noting that the points shown in Table 3 correspond to a number of material historical events. Many of these points relate to events around the 2008 global financial crisis. For example, on September 7 2008, the US Treasury Department announced to take over Fannie Mae and Freddie Mac. On 3 October 2008, the Bush administration signed a total of up to 700 billion dollars in a financial rescue plan. This shows the effectiveness of our procedure in identifying potentially influential observations to improve modeling outcomes.

5. Conclusions

In this paper, we have researched the influence diagnostics in the AR(p) model under skew-normal distributions. We have obtained the normal curvatures for the model under four perturbation schemes, including the newly proposed perturbation of skewness. We have conducted a Monte Carlo simulation study to obtain approximate benchmark values for determining those possible influential points, and use them to analyze the weekly log-returns of Chevron. The findings outlined in this paper suggest that our local influence approach in the AR(p) model effectively identifies potentially influential observations to improve the fit of the model.

Note that we have made a new application of the matrix differential calculus developed by [31] in mathematical optimization and data analysis. In particular, this paper focuses on the detection of local (rather than global) influential observations. The extension of the current study to identify global influential observations will be considered in the future.