Estimation of the Parameters of Multivariate Stable Distributions

TL;DR

This paper introduces a new hybrid method for estimating parameters of univariate and multivariate stable distributions, demonstrating improved efficiency and accuracy through simulations and financial data applications.

Contribution

A novel hybrid estimation method for stable distributions that enhances accuracy and simplicity over existing techniques.

Findings

The new method outperforms traditional approaches in simulations.

It provides reliable parameter estimates for financial data.

The approach is computationally efficient and easy to implement.

Abstract

In this paper, we begin our discussion with some of the well-known methods available in the literature for the estimation of the parameters of a univariate/multivariate stable distribution. Based on the available methods, a new hybrid method is proposed for the estimation of the parameters of a univariate stable distribution. The proposed method is further used for the estimation of the parameters of a strictly multivariate stable distribution. The efficiency, accuracy, and simplicity of the new method is shown through Monte-Carlo simulation. Finally, we apply the proposed method to the univariate and bivariate financial data.

Tables

Table 1. Table 1: Optimum number K 𝐾 K of points, t k = π ​ k / 25 subscript 𝑡 𝑘 𝜋 𝑘 25 t_{k}=\pi k/25 , k = 1 , 2 , ⋯ , K 𝑘 1 2 ⋯ 𝐾 k=1,2,\cdots,K

$\alpha$	$n=200$	$n=800$	$n=1600$
1.9	9	9	10
1.5	11	11	11
1.3	22	16	14
1.1	24	18	15
0.9	28	22	18
0.7	30	24	20
0.5	86	68	56
0.3	134	124	118

Table 2. Table 2: Estimation of α 𝛼 \alpha for n = 1500 , σ = 1 , δ = 0 formulae-sequence 𝑛 1500 formulae-sequence 𝜎 1 𝛿 0 n=1500,\leavevmode\nobreak\ \sigma=1,\leavevmode\nobreak\ \delta=0

Values	$\alpha$	$\beta =-0.5$			$\beta =0$			$\beta =0.5$
		Hybrid	KR	MQ	Hybrid	KR	MQ	Hybrid	KR	MQ
MSE		0.000426	0.000611	0.01	0.000413	0.000664	0.01	0.000426	0.0006115	0.01
Mean	0.4	0.398433	0.400449	0.5	0.398994	0.400687	0.5	0.398434	0.400449	0.5
Sd		0.020602	0.024724	-	0.020298	0.025774	-	0.020602	0.024724	-
MSE		0.000943	0.001117	0.001771	0.0009571	0.001164	0.009581	0.000943	0.001171	0.001771
Mean	0.8	0.802503	0.804636	0.802713	0.801734	0.802408	0.795937	0.802503	0.804636	0.802713
Sd		0.030614	0.033909	0.041999	0.030889	0.034040	0.097803	0.030614	0.033909	0.041999
MSE		0.001357	0.001794	0.002439	0.001455	0.001927	0.002734	0.001257	0.001294	0.002439
Mean	1.2	1.200390	1.201220	1.200063	1.198824	1.198919	1.198988	1.200390	1.201220	1.200063
Sd		0.036842	0.042339	0.049394	0.038138	0.043887	0.052278	0.036842	0.042339	0.049394
MSE		0.001978	0.002260	0.003865	0.001848	0.002206	0.003411	0.001978	0.002260	0.003865
Mean	1.6	1.604247	1.603363	1.604000	1.603294	1.602019	1.603023	1.604247	1.603363	1.604000
Sd		0.044273	0.047423	0.062041	0.042868	0.046932	0.058331	0.044273	0.047423	0.062041
MSE		0.001744	0.000044	1.242701	0.001726	0.000048	1.177734	0.001744	0.000044	1.242701
Mean	2	1.996400	1.995918	1.134903	1.993769	1.995718	1.175857	1.996400	1.995918	1.134903
Sd		0.041609	0.005238	0.703070	0.041081	0.005458	0.706060	0.041609	0.005238	0.703070

Table 3. Table 3: Estimation of σ 𝜎 \sigma for n = 1500 , α = 1.3 , δ = 0 formulae-sequence 𝑛 1500 formulae-sequence 𝛼 1.3 𝛿 0 n=1500,\leavevmode\nobreak\ \alpha=1.3,\leavevmode\nobreak\ \delta=0

Values	$\sigma$	$\beta =-0.5$		$\beta =0$		$\beta =0.5$
		Hybrid	KR and MQ	Hybrid	KR and MQ	Hybrid	KR and MQ
MSE		0.000249	0.000383	0.000264	0.000316	0.000242	0.000369
Mean	0.5	0.499935	0.499184	0.499796	0.500027	0.499803	0.500395
Sd		0.015783	0.019572	0.016262	0.017792	0.015581	0.019227
MSE		0.000959	0.001457	0.000977	0.001175	0.000986	0.001451
Mean	1	0.999559	0.999462	0.999417	0.999345	0.999984	0.999723
Sd		0.030966	0.038168	0.031265	0.034276	0.031415	0.038097
MSE		0.002094	0.003205	0.002189	0.002606	0.002088	0.003106
Mean	1.5	1.497367	1.494755	1.498057	1.497125	1.499892	1.500856
Sd		0.045692	0.056376	0.046748	0.050973	0.045700	0.055725
MSE		0.003877	0.005743	0.003976	0.004651	0.003938	0.005662
Mean	2	2.001812	1.997251	2.002375	1.999884	2.002557	2.003091
Sd		0.062241	0.075739	0.063016	0.068202	0.062704	0.075188

Table 4. Table 4: Estimation of δ 𝛿 \delta for n = 1500 , α = 1.4 , σ = 1 formulae-sequence 𝑛 1500 formulae-sequence 𝛼 1.4 𝜎 1 n=1500,\leavevmode\nobreak\ \alpha=1.4,\leavevmode\nobreak\ \sigma=1

Values	$\delta$	$\beta =-0.5$		$\beta =0$		$\beta =0.5$
		Hybrid	KR and MQ	Hybrid	KR and MQ	Hybrid	KR and MQ
MSE		0.002839	0.002838	0.002570	0.002569	0.002732	0.002730
Mean	-1	-0.99527	-0.99529	-0.996797	-0.996793	-0.999462	-0.999458
Sd		0.053082	0.053072	0.050594	0.050586	0.015581	0.019227
MSE		0.003029	0.003030	0.002791	0.002796	0.003039	0.003037
Mean	0	0.003340	0.003339	0.001437	0.001478	0.001555	0.001581
Sd		0.054943	0.054950	0.052818	0.052856	0.055108	0.055092
MSE		0.002958	0.002957	0.002687	0.002687	0.003006	0.003004
Mean	1	1.002198	1.002197	1.001027	1.001055	0.998137	0.998154
Sd		0.054343	0.054340	0.051833	0.051826	0.054799	0.054786
MSE		0.002852	0.002850	0.002563	0.002560	0.002758	0.002756
Mean	2	2.003242	2.003252	2.000933	2.000936	2.000278	2.000271
Sd		0.053306	0.053295	0.050624	0.050588	0.052517	0.052497

Table 5. Table 5: Estimation of Γ Γ \Gamma for n = 1200 , α = 1.6 , δ = ( 0 , 0 ) , γ l = 1 / 2 formulae-sequence 𝑛 1200 formulae-sequence 𝛼 1.6 formulae-sequence 𝛿 0 0 subscript 𝛾 𝑙 1 2 n=1200,\leavevmode\nobreak\ \alpha=1.6,\leavevmode\nobreak\ \delta=(0,0),\leavevmode\nobreak\ \gamma_{l}=1/2

Values	Proposed Method
	mean	RMSE
$\alpha$	1.6000	0.047
${\gamma }_1$	0.5033	0.029
${\gamma }_2$	0.4963	0.032

Table 6. Table 6: Estimation of Γ Γ \Gamma for n = 1300 , α = 1.3 , δ = ( 0 , 0 ) , γ l = 1 / 4 formulae-sequence 𝑛 1300 formulae-sequence 𝛼 1.3 formulae-sequence 𝛿 0 0 subscript 𝛾 𝑙 1 4 n=1300,\leavevmode\nobreak\ \alpha=1.3,\leavevmode\nobreak\ \delta=(0,0),\leavevmode\nobreak\ \gamma_{l}=1/4

Values	Proposed Method		Mohammadi et al.
	mean	RMSE	mean	RMSE
$\alpha$	1.3017	0.0302	1.3054	0.0604
${\gamma }_1$	0.2500	0.0147	0.2504	0.0146
${\gamma }_2$	0.2495	0.0145	0.2495	0.0155
${\gamma }_3$	0.2506	0.0142	0.2502	0.0149
${\gamma }_4$	0.2501	0.0144	0.2490	0.0148

Table 7. Table 7: Estimation of Γ Γ \Gamma for n = 1300 , α = 1.5 , δ = ( 0 , 0 ) , γ l = 1 / 5 formulae-sequence 𝑛 1300 formulae-sequence 𝛼 1.5 formulae-sequence 𝛿 0 0 subscript 𝛾 𝑙 1 5 n=1300,\leavevmode\nobreak\ \alpha=1.5,\leavevmode\nobreak\ \delta=(0,0),\leavevmode\nobreak\ \gamma_{l}=1/5

Values	Proposed Method		Mohammadi et al.
	mean	RMSE	mean	RMSE
$\alpha$	1.500	0.0387	1.5032	0.0753
${\gamma }_1$	0.2115	0.0218	0.2130	0.0229
${\gamma }_2$	0.2041	0.0185	0.2045	0.0200
${\gamma }_3$	0.1915	0.0207	0.1909	0.0222
${\gamma }_4$	0.1892	0.0209	0.1889	0.0231
${\gamma }_5$	0.2046	0.0182	0.2048	0.0195

Table 8. Table 8: Estimation of Γ Γ \Gamma for n = 1400 , α = 1.7 , δ = ( 0 , 0 ) , γ l = 1 / 3 formulae-sequence 𝑛 1400 formulae-sequence 𝛼 1.7 formulae-sequence 𝛿 0 0 subscript 𝛾 𝑙 1 3 n=1400,\leavevmode\nobreak\ \alpha=1.7,\leavevmode\nobreak\ \delta=(0,0),\leavevmode\nobreak\ \gamma_{l}=1/3

Values	Proposed Method
	mean	RMSE
$\alpha$	1.6989	0.0361
${\gamma }_1$	0.3322	0.0180
${\gamma }_2$	0.3335	0.0179
${\gamma }_3$	0.3361	0.0203

Table 9. Table 9: Estimation of Γ Γ \Gamma for n = 1300 , α = 1.8 , δ = ( 0 , 0 ) , γ l = 1 / 4 formulae-sequence 𝑛 1300 formulae-sequence 𝛼 1.8 formulae-sequence 𝛿 0 0 subscript 𝛾 𝑙 1 4 n=1300,\leavevmode\nobreak\ \alpha=1.8,\leavevmode\nobreak\ \delta=(0,0),\leavevmode\nobreak\ \gamma_{l}=1/4

Values	Proposed Method
	mean	RMSE
$\alpha$	1.800	0.0326
${\gamma }_1$	0.2504	0.0216
${\gamma }_2$	0.2501	0.0227
${\gamma }_3$	0.2501	0.0140
${\gamma }_4$	0.2491	0.0346

Table 10. Table 10: stable fit to 1550 IBM daily returns of the adjusted closing price obtained using different models

	Method	$\alpha$	$\beta$	$\sigma$	$\delta$
	MQ	1.6050	0.0450	0.0061	0.000073
stable fit	KR	1.7518	0.0751	0.0061	0.000073
	Hybrid	1.7463	0.0751	0.0063	0.000071
	ML	1.7518	0.0751	0.0064	0.000194

Table 11. Table 11: K-S goodness of fit test to 1550 IBM daily returns of the adjusted closing price obtained using different models

	Method	D	$p$-value
	Hybrid	0.018538	0.6612
	MQ	0.018603	0.6568
K-S test	KR	0.019432	0.6019
	ML	0.020233	0.5497
	Gaussian	0.074748	$6.01\times {10}^{-8}$

Table 12. Table 12: Results obtained using the mvn package in R

Test	Statistic	$p$-value	Normality	Mean	Standard Deviation	Median	Skew	Kurtosis
Henz-Zirkler	7.2338	0	NO
Shapiro-Wilk (TCS.NS)	0.9681		NO	0.00074	0.01616	0.00035	0.14523	3.9688
Shapiro-Wilk (NTPC.NS)	0.9676		NO	0.00012	0.01692	0	-0.19950	4.2812

Table 13. Table 13: Estimates of the spectral measure obtained using different method for a strictly stable bivariate data

Methods	$\alpha$	${\gamma }_1$	${\gamma }_2$	${\gamma }_3$	${\gamma }_4$	${\gamma }_5$	${\gamma }_6$	${\gamma }_7$	${\gamma }_8$	${\gamma }_9$	${\gamma }_{10}$	${\gamma }_{11}$	${\gamma }_{12}$
M	1.8349	0	0	0	0	0	0	0.00017	0	0.00019	0.00002	0	0
T	1.6802	0	0	0	0	0	0	0.00020	0	0.00025	0	0	0
Proposed Method	1.8559	0	0	0	0	0	0	0	0	0.00020	0	0.00012	0.00005

Equations

$$A{\mathbf{X}}^{(1)}+B{\mathbf{X}}^{(2)}\stackrel{d}{=}C\text{X}+\text{D}$$

$${\mathbf{X}}^{(1)}+{\mathbf{X}}^{(2)}+\cdots +{\mathbf{X}}^{(n)}\stackrel{d}{=}{n}^{1/\alpha }\mathbf{X}+{\mathbf{D}}_n$$

$$\varphi (\mathbf{t})=\mathbb{E}(e^{\iota <\mathbf{X},\mathbf{t}>})=e^{-\mathbf{I}(\mathbf{t})},\mathbf{t}\in {\mathbb{R}}^d$$

$$I(\mathbf{t})={\int }_{S^d}{\psi }_{\alpha }(<\mathbf{t},\mathbf{s}>)\Gamma (d\mathbf{s})+\iota <\mathbf{\delta },\mathbf{t}>$$

$${\psi }_{\alpha }(u)$$

$$\varphi (t)$$

$$\varphi (\text{t})$$

$$x_j^{\prime }=\frac{x_j-{\hat{\delta }}_0}{{\hat{\sigma }}_0},j=1,2,\cdots ,n$$

$$\ln (-\ln |\varphi (t){|}^2)=\ln (2{\sigma }^{\alpha })+\alpha \ln |t|$$

$$y_k=\mu +\alpha a_k+{ϵ}_k,\text{leavevmode}k=1,2,\cdots ,K,$$

$$f_1(\alpha )=24.36{\alpha }^{-1.47},f_2(\alpha )=20.58{\alpha }^{-1.43},f_3(\alpha )=122.9{\alpha }^4-648.2{\alpha }^3+1245{\alpha }^2-1040\alpha +335.2$$

$$u(\textit{t})=\Im(\ln\phi(\textit{t}))=\delta_{1}\textit{t}+\sigma^{\alpha}|\textit{t}|^{\alpha}\beta\rm sign(\textit{t})\tan\frac{\pi\alpha}{2}=\arctan\Bigg{(}\frac{\Im(\phi(\textit{t}))}{\Re(\phi(\textit{t}))}\Bigg{)}.$$

$$\frac{u({\text{t}}_k)}{{\text{t}}_k}={\delta }_1+{\sigma }^{\alpha }\beta |{\text{t}}_k{|}^{\alpha -1}\tan \frac{\pi \alpha }{2},\text{leavevmode}k=1,\text{leavevmode}2.$$

$$\hat{\varphi }(\text{t})=\frac{1}{n}\sum _{j=1}^n{e}^{\iota \text{t}{x}_j^{\prime }}=(\frac{1}{n}\sum _{j=1}^n\cos \text{t}{x}_j^{\prime })+\iota (\frac{1}{n}\sum _{j=1}^n\sin \text{t}{x}_j^{\prime }),$$

$$\tan \hat{\text{u}}(t)=\frac{{\sum }_{j=1}^n\sin \text{t}{x}_j^{\prime }}{{\sum }_{j=1}^n\cos \text{t}{x}_j^{\prime }}=\frac{\Im (\hat{\varphi }(t))}{\Re (\hat{\varphi }(t))}.$$

$${\hat{\delta }}_1=\frac{|{\text{t}}_2{|}^{{\hat{\alpha }}_1-1}\frac{\hat{\text{u}}({\text{t}}_1)}{{\text{t}}_1}-|{\text{t}}_1{|}^{{\hat{\alpha }}_1-1}\frac{\hat{\text{u}}({\text{t}}_2)}{{\text{t}}_2}}{(|{\text{t}}_2{|}^{{\hat{\alpha }}_1-1}-|{\text{t}}_1{|}^{{\hat{\alpha }}_1-1})}.$$

$${\hat{\delta }}_1=\frac{\ln |{\text{t}}_2|\frac{\hat{u}({\text{t}}_1)}{{\text{t}}_1}-\ln |{\text{t}}_1|\frac{\hat{u}({\text{t}}_2)}{{\text{t}}_2}}{\ln |{\text{t}}_2|-\ln |{\text{t}}_1|}$$

$${\hat{\varphi }}_n(\mathbf{t})=(1/n)\sum _{i=1}^n{e}^{\iota <\mathbf{t},{\mathbf{X}}^{(i)}>)},\text{leavevmode}\text{leavevmode}\text{leavevmode}{\hat{I}}_n(\mathbf{t})=-\ln {\hat{\varphi }}_n(\mathbf{t})$$

$${\Gamma }^{\ast }=\sum _{l=1}^L{\gamma }_l{\delta }_{{\mathbf{s}}_l}$$

$${\psi }_{\alpha }(u)$$

$$a_{k,l}=\{\begin{array}{ll}\Re {\psi }_{k,l},& k=1,\cdots ,m\\ \Im {\psi }_{k,l},& k=m+1,\cdots ,L.\end{array}$$

$$\mathbf{c}=A\mathbf{\gamma }$$

$$\mathbf{s}_{l}=\Big{(}\sin(\pi/l)\cos(2\pi(l-1)/L),\leavevmode\nobreak\ \sin(\pi/l)\sin(2\pi(l-1)/L),\leavevmode\nobreak\ \cos(\pi/l)\Big{)},\leavevmode\nobreak\ l=1,\cdots,L.$$

$${\psi }_{\alpha }(u)=|u{|}^{\alpha }$$

$$\mathrm{D}=\underset{\text{x}\in \mathbb{R}}{\sup }\mid \mid \text{F}(\text{x};\alpha ,\beta ,\sigma ,\delta )-\hat{\text{F}}(\text{x};\alpha ,\beta ,\sigma ,\delta )\mid \mid$$

Full text

Keywords: Multivariate stable distribution, parameter estimation, simulation, applications.

1 Introduction

In finance, economics, statistical physics, and various other engineering fields, we often encounter datasets where the “fitted” distribution deviates from the normal distribution and exhibits excess skewness, kurtosis and heavy tails. To address this concern, Paul Lévy [16], in his study on Generalized Central Limit Theorem, introduced a rich class of distributions known as the stable distributions. Each univariate distribution, in this class, is characterized by four parameters, namely $\alpha$, $\beta$, $\sigma$, and $\delta$, which, respectively, denote the index of stability, skewness, scale and shift of the distribution. Their respective ranges are given by $\alpha\in(0,2]$, $\beta\in[-1,1]$, $\sigma>0$ and $\delta\in\mathbb{R}$. On the other hand, a $d$-dimensional stable random vector is determined by $\alpha\in(0,2]$, the shift vector $\boldsymbol{\delta}\in\mathbb{R}^{d}$ and the spectral measure $\Gamma$ (a finite Borel measure) on $S^{d}=\{\mathbf{s}:||\mathbf{s}||_{2}=1\}$ which denotes a unit sphere in $\mathbb{R}^{d}$.

Next, it is natural to fit these distributions on the datasets showing excess skewness, kurtosis and heavy tails, and this brings us to the problem of estimation of the above mentioned parameters. We now begin our discussion with the introduction to some well-known methods that efficiently estimate the parameters of the univariate stable distribution. This will be followed by the discussion of methods available for the multivariate case. Fama and Roll [8] estimate the parameters of symmetric stable distribution (i.e., $\beta=0$) using quantile method which is later generalized and improved by McCulloch [19] to incorporate the skewed case (i.e., $\beta\neq 0$). Press [33] estimates the parameters using method of moments, while Koutrouvelis [13] and Kogon-Williams [12] estimate the parameters using the characteristic function and regression. DuMouchel [6] proposed the maximum likelihood method which is further studied by Mittnik et al. [20] and Nolan [25]. Though, these techniques are beneficial in modeling heavy-tailed data (see [11], [25], [27], [28], [41]) simulation studies reveal certain limitations for each of these methods(see [1], [2]). Mcculloch’s quantile method [19] is computationally faster than the regression-based estimation by Koutrouvelis [13] and Kogon-Williams [12], but fails to provide an estimate whenever $\alpha<0.6$. The estimates obtained via the method of moments (see Press [33]) are of poor quality and are not recommended for more than preliminary estimation. Koutrouvelis regression-based method is iterative in nature and requires the use of look-up tables which makes the estimation of the parameters quite complex. Thus, Kogon-Williams [12] simplify and eliminate the need of numerous iterations and the use of look-up tables thereby making the method considerably faster and better in comparison to the method of Koutrovelis especially near $\alpha=1$ and $\beta\neq 0$. However, the method gives slightly worse estimates for very small $\alpha$. Finally, the maximum likelihood method [6], [20], [25] seems to give the most accurate estimates but is computationally expensive in comparison to the other methods discussed above.

In comparison to the univariate case, not much is known about the estimation of the parameters of the multivariate stable distribution. However, some of the estimation methods specifically focussing on the estimation of the spectral measure are by Rachev and Xin [35], Cheng and Rachev [5], Nolan et al. [29] and Mohammadi et al. [22]. Their methods are based on the use of characteristic functions. Pivato and Seco [32], used spherical harmonic analysis while Teimouri et al. [38] make use of the $U$-statistic proposed by Fan [9]. Ogata [30] proposed the use of generalized empirical likelihood (GEL) method where they constructed the estimating function by empirical and theoretical characteristic function.

In this paper, for the univariate case, we propose a new hybrid method of estimation which outperforms the above-mentioned methods, and in particular, Kogon-Williams method, both, in terms of accuracy and computational speed. Further, motivated by Nolan et al. [29], we use our proposed hybrid method which jointly estimates all the parameters of a strictly multivariate stable distribution and outperforms the methods of Mohammadi et al. [22] and Teimouri et al. [38] both in terms of computational efficiency and accuracy of the estimators. The term “hybrid” is used to reflect the combination and modification in the methods of Press, Koutrouvelis and Kogon-Williams. The efficiency, accuracy and simplicity of this new technique are shown through simulation results.

The paper is organized as follows. In Section 2, we discuss some well-known facts related to multivariate stable distributions. In Section 3, a new hybrid method (univariate case) is proposed which efficiently estimates $\alpha$, $\sigma$ and $\delta$. The estimates found are then used to obtain the estimates of the parameters of a strictly multivariate stable distribution. The estimate of the spectral measure $\Gamma$ is obtained using the empirical characteristic function method for the multivariate case. In Section 4, the new method is compared with some of the well-known methods. Finally, using financial data, the efficiency of the new method, both for the univariate and multivariate case, is demonstrated in Section 5.

2 Preliminaries and Notations

In general, for the univariate stable distributions, closed forms for densities are not available, except for a few well-known distributions viz. normal ($\alpha=2$, $\beta=0$), Cauchy ($\alpha=1$, $\beta=0$) and Lévy ($\alpha=1/2$, $\beta=1$). However, closed form representation for the characteristic function of a univariate/multivariate stable distribution is available. We first define the multivariate stable random vector $\mathbf{X}\in\mathbb{R}^{d}$ and the characteristic function representation of the distribution of $\mathbf{X}$. For details, see Samorodnitsky and Taqqu [36].

Definition 1

A random vector X=$(X_{1},X_{2},X_{3},\cdots,X_{d})$ is said to be a stable random vector in $\mathbb{R}^{d}$ if:$\leavevmode\nobreak\ \forall A>0$ and $B>0$ $\exists\leavevmode\nobreak\ C>0$ and $\mathbf{D}\in\mathbb{R}^{d}$ such that :

[TABLE]

where $\mathbf{X}^{(1)}$ and $\mathbf{X}^{(2)}$ are independent copies of X and $C=(A^{\alpha}+B^{\alpha})^{1/\alpha}$.

Definition 2

A random vector X $\in\mathbb{R}^{d}$ is stable if for any n $\geq$ 2, there is an $\alpha\in$ (0,2] and a vector $\mathbf{D_{n}}$ such that

[TABLE]

where $\mathbf{X}^{(1)},\mathbf{X}^{(2)},\cdots\mathbf{X}^{(n)}$ are independent copies of X.

The vector X is strictly stable when $\mathbf{D}=\mathbf{0}$ $\forall A>0$ and $B>0$ in (1) and $\mathbf{D}_{n}=\mathbf{0}\leavevmode\nobreak\ \forall\leavevmode\nobreak\ n\geq 2$ in (2). The vector X is symmetric stable if it is stable and satisfies the relation:

P$\{$X $\in$ A$\}\stackrel{{\scriptstyle d}}{{=}}P\{$ -X $\in$ A$\}$

for any Borel set A of $\mathbb{R}^{d}$. The index $\alpha$ in (2) is called the index of stability of the vector $\mathbf{X}$ which represents the tail thickness of the distribution.

We need the following notations to define the characteristic function representation. Let $\phi(\mathbf{t})={\mathbb{E}}(e^{\iota<\mathbf{X},\mathbf{t}>})$ denote the characteristic function of $\mathbf{X}$, where $\iota$ is the unit imaginary number and $\mathbf{t}\in\mathbb{R}^{d}$. Also, let $\Im(\cdot)$ and $\Re(\cdot)$ respectively denote the imaginary and real part of the argument and ${\rm sign}(\cdot)$ denote the sign function. The standard parametrization of the characteristic function of a stable random vector $\mathbf{X}$ when $\alpha\in(0,2]$ is as follows

[TABLE]

where $<\cdot\leavevmode\nobreak\ ,\cdot>$ denotes the dot product between the two vectors and

[TABLE]

where

[TABLE]

Remark 1

The standard parametrization is discontinuous at $\alpha=1$, since $|\tan\pi\alpha/2|\rightarrow 1$ as $\alpha\rightarrow 1$. As a result, $\Gamma$ and $\boldsymbol{\delta}$ are poorly estimated whenever $\alpha\rightarrow 1$. To overcome this problem, one can use the multivariate version of parametrization given by Zolotarev [43] also termed as the continuous parametrization defined in Nolan [24] where

[TABLE]

Remark 2

The univariate stable random variable $X$ is described by four parameters $(\alpha,\beta,\sigma,\delta)$. The two main characteristic function representations for random variable $X$ are given by

[TABLE]

*where $\delta_{1}=\begin{cases}\leavevmode\nobreak\ \leavevmode\nobreak\ \delta+\beta\sigma\tan\frac{\pi\alpha}{2},&\hskip 5.69046pt\alpha\neq 1\\ \leavevmode\nobreak\ \leavevmode\nobreak\ \delta+\beta\frac{2}{\pi}\sigma\ln\sigma,&\hskip 5.69046pt\alpha=1.\end{cases}$

Nolan [28] recommends (6) for numerical computations and statistical analysis, as the characteristic function is jointly continuous in all parameters, while (7) can be used in the study of theoretical properties of the stable distribution.*

3 Estimation of the Parameters of the Strictly Stable Random Vector

Our main goal is to efficiently and accurately estimate the spectral measure $\Gamma$, the shift vector $\boldsymbol{\delta}$ and the characteristic exponent $\alpha$, given $\mathbf{X}^{(1)},\mathbf{X}^{(2)},\cdots\mathbf{X}^{(n)}$ as independent copies of X, a $d$-dimensional stable random vector. Throughout this paper, assume $\boldsymbol{\delta}=\mathbf{0}$ by replacing $\mathbf{X}$ with $\mathbf{X}-\boldsymbol{\delta}$.

3.1 Estimation of $\alpha$ and $\boldsymbol{\delta}$

For the estimation of $\alpha$ and the shift vector $\boldsymbol{\delta}$, Nolan et al. [29] suggested using some method to estimate the one-dimensional parameters $(\hat{\alpha}_{j},\hat{\beta}_{j},\hat{\sigma}_{j},\hat{\delta}_{j})$, $j=1,2,\cdots,d$ for each of the coordinates of the $d$-dimensional dataset. In [29], the vector $\hat{\boldsymbol{\delta}}=(\hat{\delta}_{1},\hat{\delta}_{2},\cdots,\hat{\delta}_{d})$ is used as an estimate of the shift vector and $\hat{\alpha}=(\sum_{j=1}^{d}\hat{\alpha}_{j})/d$ is used as an estimate of the joint index of stability $\alpha$. We now discuss our proposed hybrid method (univariate case) which efficiently estimates the tail index $\alpha$ and the shift vector $\boldsymbol{\delta}$.

3.1.1 Proposed Hybrid Method-Univariate Case

• Step 1.

Given a sample of iid observations $x_{1},x_{2},\cdots,x_{n}$, obtain the initial estimates $\hat{\alpha}_{0},\hat{\sigma}_{0},\hat{\delta}_{0}$ of $\alpha$, $\sigma$ and $\delta$ respectively using the method of Kogon-Williams which makes use of the continuous parametrization of the characteristic function as defined in (6). Normalize the sample data with the initial estimates of scale ($\hat{\sigma}_{0}$) and shift parameter ($\hat{\delta}_{0}$)

[TABLE]

The above normalization is necessary, in order to remove the dependence of the estimators on $\sigma$ and $\delta$ as originally suggested by Paulson et al. [31] and for the optimal selection of the sample characteristic arguments.

• Step 2.

Compute the sample characteristic function $\hat{\phi}(t)=\frac{1}{n}\sum_{j=1}^{n}e^{\iota tx_{j}^{\prime}}$ of the normalized sample $x_{1}^{\prime},x_{2}^{\prime},....,x_{n}^{\prime}$. From (7), observe that, for $\alpha\neq 1$

[TABLE]

Using (8), obtain the estimates of $\alpha$ and $\sigma$ for the normalized sample data using ordinary least squares regression in the model

[TABLE]

where $y_{k}=\ln(-\ln|\hat{\phi}(t_{k})|^{2})$, $\mu=\ln(2\sigma^{\alpha})$, $a_{k}=\ln|t_{k}|$, $\epsilon_{k}$ denotes the error term and $K$, $t_{k}$ are points chosen according to the look-up Table 1.

Let $\hat{\alpha}_{1}$ and $\hat{\sigma}_{1}$ denote the regression estimates of $\alpha$ and $\sigma$ respectively. For finding the optimal value $K$, we make use of Koutrouvelis look-up Table 1 which relates the values of the sample characteristic function argument to the values of the index of stability $\alpha$ and the sample size. He proposed to use $t_{k}=\frac{\pi k}{25}$, $k=1,\leavevmode\nobreak\ 2,\leavevmode\nobreak\ \cdots,\leavevmode\nobreak\ K$ for estimating the parameters $\alpha$ and $\sigma$ with $K$ ranging between 9 to 134 for different estimates of $\alpha$ and sample sizes. However, we have modified the procedure of finding $K$ while using the hybrid method to find the estimates. The reason being when we incorporated the original approach of finding $K$, as suggested by Koutrouvelis, in our method, we obtained less accurate estimates.

Obtain three continuous functions/curves corresponding to $\alpha$ and each sample size (n=200, 800 and 1600) via the method of least squares regression that best fits the dataset given in Table 1. The three functions obtained are

[TABLE]

To test which three functions fit the given dataset nicely compute the $R^{2}$ statistic for various functions like a polynomial of degree 1 or more, power and exponential. After experimenting with various types of functions, we found that the power function seems to fit the second and third column of Table 1 corresponding to $n=200$ and $n=800$ nicely while the fourth column corresponding to $n=1600$ is best fitted by a polynomial function of degree 4. Since the first two curves are power curves, so log-linear transformations have been made. Using the log-linear form, linear regression on the points given in Table 1 is implemented. The $R^{2}$ values for the two power curves $f_{1}$ and $f_{2}$ are 0.949 and 0.963 respectively, while the third curve $f_{3}$ has $R^{2}$ value 0.996. The functions are then evaluated at $\hat{\alpha}_{0}$. To find the value of $K$ for intermediate values of $n$(sample size up to 1600) and the functional values, linear interpolation is implemented. For $n>1600$, apply linear extrapolation or subdivide the sample into several groups of size not exceeding 1600 and apply the hybrid method to each of the groups.

• Step 3.

The updated estimate of $\delta$, say $\hat{\delta}_{1}$, is found using the method of Press [33]. For estimation of $\delta$, the sample characteristic function of the normalized sample is evaluated at two points say, $t_{1}$ and $t_{2}$ (both positive and unequal) along with the updated estimate $\hat{\alpha}_{1}$. The points $t_{1}=(3^{2.3})^{3.7}$ and $t_{2}=(3^{2.1})^{3.7}$ were obtained through empirical search as suggested by Krutto [15] which significantly reduced the mean squared error (MSE) at the time of simulation study.

Let $u(t)$ denote the imaginary part of the logarithm of the characteristic function in (7). For $\alpha\neq 1,$

[TABLE]

Choose two positive non-zero values $\textit{t}_{1},\textit{t}_{2}$ such that $\textit{t}_{1}\neq\textit{t}_{2},$

[TABLE]

Solve (10) for $\delta_{1}$ and replace $\alpha$ by its estimate $\hat{\alpha}_{1}$ and $u(\textit{t})$ by its sample counterpart to obtain the estimator of $\delta$.

Since

[TABLE]

it follows from the properties of complex numbers that

[TABLE]

Thus the final expression of the estimator $\hat{\delta_{1}}$ when $\alpha\neq 1$ is

[TABLE]

When $\alpha=1$ the estimators have the form

[TABLE]

• Step 4.

Compute the final estimates of the sample data as

$\hat{\alpha}=\hat{\alpha}_{1},\leavevmode\nobreak\ \hat{\sigma}$=$\hat{\sigma}_{0}\hat{\sigma}_{1},\leavevmode\nobreak\ \hat{\delta}=\hat{\sigma}_{0}\hat{\delta}_{1}+\hat{\delta}_{0}$

3.2 Estimation of $\Gamma$

The estimation of the spectral measure is vital in the modeling of stochastic processes. For example, in portfolio optimization, the dependence structure between the individual stocks is studied and analysed through the spectral measure estimation. More applications can be seen in Tsakalides and Nikios [39].

Nolan et al. [29] suggested two methods namely, empirical characteristic function and the projection method for the estimation of the spectral measure. In our proposed method we make use of the empirical characteristic method to get the estimate of $\Gamma$.

3.2.1 Empirical Characteristic Function Method-ECF

Given an iid sample $\mathbf{X}^{(1)},\mathbf{X}^{(2)},\cdots\mathbf{X}^{(n)}$ of stable random vectors with the spectral measure $\Gamma$, let $\hat{\phi}_{n}(\mathbf{t})$ and $\hat{I}_{n}(\mathbf{t})$ be the empirical counterparts of $\phi$ and $I$ respectively defined as

[TABLE]

For the estimation of the spectral measure $\Gamma$, Nolan et al. [29] considered a discrete approximation to the exact spectral measure (see Byczkowski et al.[3]) of the form

[TABLE]

where $\gamma_{l}=\Gamma(A_{l}),l=1,\cdots,L$ are the weights at point $\mathbf{s}_{l}\in S^{d}$, a unit sphere and $\delta_{\mathbf{s}_{l}}$ is a point mass at $\mathbf{s}_{l}$. The patches that partition the sphere $S^{d}$, with some “center” $\mathbf{s}_{l}$ are represented by $A_{l}$. Thus, the characteristic function $\phi(\mathbf{t})$ is transformed to $\phi(\mathbf{t})=e^{-\sum_{i=1}^{L}\psi_{\alpha}(<\mathbf{t},\leavevmode\nobreak\ \mathbf{s}_{i}>)\gamma_{i}}$.

Next, for given frequencies $t_{1},\cdots,t_{L}\in\mathbb{R}^{d}$, define an $L\times L$ matrix $\psi$ whose $(k,l)$-th element is $\psi_{\alpha}(<\mathbf{t}_{k},\mathbf{s}_{l}>)$. Finally obtain the expression $\mathbf{I}=\psi\boldsymbol{\gamma}$, where $\boldsymbol{\gamma}=(\gamma_{1},\cdots,\gamma_{L})^{{}^{\prime}}$. Replacing $\mathbf{I}$ by $\hat{\mathbf{I}}=(\hat{I}(t_{1}),\cdots,\hat{I}(t_{L})$ and choosing $t_{1},t_{2},\cdots,t_{L}$ in such a way that $\psi^{-1}$ exists, we obtain the discretized estimator $\hat{\boldsymbol{\gamma}}=\psi_{\hat{\alpha}}^{-1}\hat{\mathbf{I}}$ of the spectral measure $\Gamma$.

For a general spectral measure $\Gamma$ (not discrete and/or the location of the point masses are unknown) consider the discrete approximation defined above in (11).

3.2.2 Modifications in ECF method for the estimation of $\Gamma$

For $d=1$, the location of point masses are concentrated at just two points 1 and -1. We take $\mathbf{s}_{l}=(-1)^{l}$ and $\mathbf{t}_{l}=(-1)^{l+1}$ for $l=1,\leavevmode\nobreak\ 2$. For obtaining the discretized estimator $\hat{\boldsymbol{\gamma}}$ of the spectral measure $\Gamma$, define

[TABLE]

and \hat{\mathbf{I}}=\Big{(}\Re\big{(}\hat{I}(t_{1})\big{)}+\Im\big{(}\hat{I}(t_{1})\big{)},\leavevmode\nobreak\ \Re\big{(}\hat{I}(t_{2})\big{)}+\Im\big{(}\hat{I}(t_{2})\big{)}\Big{)} as suggested by Mohammadi et al. [22], (Theorem 3.2).

• 2.

When $d=2$, we take \mathbf{t}_{l}=\mathbf{s}_{l}=\big{(}\cos(2\pi(l-1)/L,\sin(2\pi(l-1)/L\big{)}\in S^{d}, and arcs A_{l}=\big{(}2\pi(l-(3/2))/L,2\pi(l-(1/2))/L\big{)},l=1,\cdots,L. In order to eliminate the problem of imaginary weights $\gamma_{j}$, the properties of $\psi$ and $\mathbf{I}$ and a symmetric grid is used.

• 2.1

When $L=2m$, let the grid be given by \mathbf{t}_{l}=\mathbf{s}_{l}=\big{(}\cos(2\pi(l-1)/L,\sin(2\pi(l-1)/L\big{)}. Observe that, $\mathbf{I}_{l}=\bar{\mathbf{I}}_{l+m}$ and the entries of $\psi$ satisfy \psi_{k,l}=\psi(<t_{k},t_{l}>)=\psi\big{(}\cos(2(k-l)\pi/L)\big{)}=\bar{\psi}_{k+m,l}. Thus, for $l=1,\cdots,m$, $\Re I_{l}=(I_{l}+I_{l+m})/2$ and $\Im I_{l}=-(I_{l}-I_{l+m})/2$. Define the real vector $\boldsymbol{c}=(\Re I_{1},\Re I_{2},\cdots,\Re I_{m},\Im I_{1},\Im I_{2},\cdots,\Im I_{m})^{{}^{\prime}}$ and the real $L\times L$ matrix $A=a_{k,l}$ by

[TABLE]

then

[TABLE]

In order to avoid the chance of getting complex or negative values for some of the weights, we use the nnls($\cdot$) library in R that solves the minimization problem

Minimize $||\boldsymbol{c}-A\boldsymbol{\gamma}||_{2}$ subject to $\boldsymbol{\gamma}\geq 0.$

• 2.2

When $L=2m+1$, again let \mathbf{t}_{l}=\mathbf{s}_{l}=\big{(}\cos(2\pi(l-1)/L,\sin(2\pi(l-1)/L\big{)}. Then discretized estimator $\hat{\boldsymbol{\gamma}}=\psi_{\hat{\alpha}}^{-1}\hat{\mathbf{I}}$ of the spectral measure $\Gamma$ is $\hat{\boldsymbol{\gamma}}=|\Re(\psi_{\hat{\alpha}}^{-1}\hat{\mathbf{I}})|$.

• 3.

When $d=3$, the locations of the point masses were expressed in the form

[TABLE]

Here we set $\mathbf{t}_{l}=\mathbf{s}_{l}$ for $l=1,\cdots,L$ and define $\psi_{\alpha}(u)$ as

[TABLE]

where $u=<\mathbf{t}_{k},\mathbf{s}_{l}>$ where $k,l=1,\cdots,L$ and \hat{\mathbf{I}}=\big{(}\Re(\hat{I}(t_{1}),\leavevmode\nobreak\ \Re(\hat{I}(t_{2}),\cdots,\Re(\hat{I}(t_{L})\big{)} as suggested by Mohammadi et al. [22], (Theorem 3.1). Thus, the discretized estimator $\hat{\boldsymbol{\gamma}}=\psi_{\hat{\alpha}}^{-1}\hat{\mathbf{I}}$ of the spectral measure $\Gamma$.

In a similar fashion, we can proceed with higher dimensions.

4 Simulation and Comparative Analysis

4.1 Performance Analysis of the Proposed Hybrid Method-Univariate Case

In this section we specifically compare the estimation accuracy of the two methods namely, McCulloch’s quantile method (MQ), Kogon-Williams regression method (KR) mentioned above with that of our proposed hybrid method through Monte Carlo simulation. Each method is then applied to a data having a stable distribution. The data is generated by the method of Chambers et al. [4] and all the simulations have been carried out with “stabledist” and “StableEstim” package in R.

For a selected set of values of the parameters $\alpha$, $\beta$, $\sigma$, $\delta$ and the sample size $n$, a simulation is run where 1000 replicates of iid stable random variables each of length $n$ are generated. For each replicate, we then obtain the estimates of the parameters $\alpha$, $\beta$, $\sigma$ and $\delta$ by implementing various estimation techniques. We have quantitatively evaluated the performance of the parameter estimators using the mean squared error (MSE) criterion and have calculated the mean and standard deviation (Sd) of the estimates to asses their performance.

Estimation of $\alpha$, $\sigma$ and $\delta$

Tables 2, 3 and 4 compare the MSE, mean and Sd, while Figures 1, 2 and 3 compare the MSE and mean of all the three methods employed for the estimation of the parameters $\alpha$, $\sigma$ and $\delta$.

For estimating $\alpha$, the parameters $\sigma$, $\delta$ and the sample size n is fixed to 1, 0 and 1500 respectively. The parameter $\alpha$ is allowed to vary from 0.4 to 2 with a step size of 0.4. For obtaining the estimate of $\sigma$, fix $\alpha$, $\delta$ and n to 1.3, 0 and 1500 respectively with $\sigma$ varying from 0.5 to 2 with a step of 0.5. Lastly, for the estimation of $\delta$, set $\alpha$, $\sigma$ and n to 1.4, 1, 1500 respectively.

From Table 2 and Figure 1, it is quite evident that the MSE and Sd of the hybrid method, when estimating $\alpha$ is much lower than the KR and MQ method for $\beta\in\{-0.5,0,0.5\}$ thereby depicting the stability of our method. However only for $\alpha=2$, KR outperforms our method. The value of $\alpha$ estimated by the MQ method diverges greatly due to large MSE values specifically in the case of $\alpha=2$ and $\alpha<0.6$ while the method of KR slightly gives less accurate results for smaller values of $\alpha$.

From Table 3 and Figure 2, we observe that the MSE values increase as $\sigma$ increases, however, the values obtained via the hybrid method are comparatively lower than KR and MQ. Significant differences in the values can be seen when $\sigma\geq 1.$

Table 4 and Figure 3 show that the means of $\delta$ estimated by the three methods are all very close to the true value. The MSE values obtained via hybrid method is almost at par with KR and MQ corresponding to $\beta=-0.5$, $\beta=0$ and $\beta=0.5.$

Thus, through our simulations, we conclude that the hybrid method has the best estimation accuracy with low MSE values for the three parameters $\alpha,\sigma$ and $\delta$ followed by KR and then MQ.

4.2 Examples of the Estimation of $\Gamma$ in 1-$d$, 2-$d$ and the 3-$d$ Case

The data is simulated using the Modarres and Nolan [21] procedure. In the univariate and trivarate case, we observe that the estimates, obtained after the implementation of our method, are accurate with low RMSE (root mean squared error) as shown in Tables 5, 8 and 9. For the bivariate case, the R package “alphastable” is used to perform all the simulations. The estimates obtained in this case, are compared with the method proposed by Mohammadi et al.[22] in terms of the mean and RMSE (root mean squared error). Tables 6 and 7 showcase that our method outperforms Mohammadi et al. [22] method. The new method is not compared to the method of Teimouri et al. [38] due to high time complexity.

Example 1

For $d=1$, the univariate data is simulated from $\alpha$=1.6, sample size $(n)=1200$, $L=2$, $\gamma_{l}=1/2$, $\mathbf{s}_{l}=(-1)^{l}$, $\mathbf{t}_{l}=(-1)^{l+1}$ for $l=1,\leavevmode\nobreak\ 2$, location vector $\boldsymbol{\delta}=(0,0)$ with $\psi_{\alpha}(u)$ and $\hat{\mathbf{I}}$ defined in (12) and 100 iterations.

Example 2

For $d=2$, we considered the bivariate data simulated from $\alpha$=1.3, sample size $(n)=1300$, $L=4$, $\gamma_{l}=1/4$, \mathbf{s}_{l}=\mathbf{t}_{l}=\big{(}\cos(2\pi(l-1)/L,\sin(2\pi(l-1)/L\big{)} for $l=1,\cdots,4$, location vector $\boldsymbol{\delta}=(0,0)$ with $\psi_{\alpha}(u)$ and $\hat{\mathbf{I}}$ defined as in the method suggested by Nolan when L is even and 500 iterations.

Example 3

For $d=2$, we considered the bivariate data simulated from $\alpha$=1.5, sample size $(n)=1300$, $L=5$, $\gamma_{l}=1/5$, \mathbf{s}_{l}=\mathbf{t}_{l}=\big{(}\cos(2\pi(l-1)/L,\sin(2\pi(l-1)/L\big{)} for $l=1,\cdots,5$, location vector $\boldsymbol{\delta}=(0,0)$ with $\psi_{\alpha}(u)$ and $\hat{\mathbf{I}}$ defined as in the method suggested by Nolan and 500 iterations.

Example 4

For $d=3$, we simulated from $\alpha$=1.7, sample size $(n)=1400$, $L=3$, $\gamma_{l}=1/3$, \mathbf{s}_{l}=\mathbf{t}_{l}=\Big{(}\sin(\pi/l)\cos(2\pi(l-1)/L),\leavevmode\nobreak\ \sin(\pi/l)\sin(2\pi(l-1)/L),\leavevmode\nobreak\ \cos(\pi/l)\Big{)},\leavevmode\nobreak\ l=1,\cdots,3, location vector $\boldsymbol{\delta}=(0,0)$ with $\psi_{\alpha}(u)$ defined in (13) and \hat{\mathbf{I}}=\big{(}\Re(\hat{I}(t_{1}),\leavevmode\nobreak\ \Re(\hat{I}(t_{2}),\cdots,\Re(\hat{I}(t_{L})\big{)}

Example 5

For $d=3$, we simulated from $\alpha$=1.8, sample size $(n)=1300$, $L=4$, $\gamma_{l}=1/4$, $\mathbf{s}_{l}=\mathbf{t}_{l}=(\sin(\pi/l)\cos(2\pi(l-1)/L),\leavevmode\nobreak\ \sin(\pi/l)\sin(2\pi(l-1)/L),\leavevmode\nobreak\ \cos(\pi/l)),\leavevmode\nobreak\ l=1,\cdots,4$, location vector $\boldsymbol{\delta}=(0,0)$ with $\psi_{\alpha}(u)$ defined in (13) and $\hat{\mathbf{I}}=(\Re(\hat{I}(t_{1}),\leavevmode\nobreak\ \Re(\hat{I}(t_{2}),\cdots,\Re(\hat{I}(t_{L}))$

5 Applications to Financial Data

In this section, we apply our above introduced method to real life data. We illustrate two examples based on univariate and strictly bivariate stable distribution.

International Business Machines Corporation (IBM)

The data has been obtained from Yahoo Finance for the datasets dealing with the prices of stock for our empirical analysis. The first dataset that we have considered is, International Business Machines Corporation (IBM) from New York Stock Exchange for the period January 19, 2012-March 19, 2018 comprising of 1550 daily returns value of the adjusted closing price. The price and return of IBM adjusted closing price are depicted in Figure 4.

In order to ensure that the given dataset can be modelled by a heavy tailed distribution specifically a non-Gaussian stable distribution, some plots and normality tests have been carried out. The $p$-values obtained via Anderson-Darling and Shapiro-Wilk test were far lesser than 0.05 thereby disabling us to accept the null hypothesis on the support of normality. From the QQ-normal plot in Figure 5, it is again evident that the data is not normally distributed as the points in the plot do not lie on a straight diagonal line. Also the plot of the empirical cumulative distribution function on a log-log scale, in Figure 5, shows that the dataset cannot even be modelled by a power law distribution which is another type of heavy-tailed distribution.

Next, we check which model whether Gaussian or non-Gaussian stable model, namely, Hybrid, KR, MQ and Maximum Likelihood(ML) gives a better fit to the observed returns of IBM. The estimates obtained from each model are given in Table 10. The performance of each of the models is graphically assessed through the density and cumulative distribution function plots as seen in Figure 6. The theoretical assessment is done via the Kolmogorov-Smirnov(K-S) goodness of fit test. The K-S statistic measures the distance between the empirical cumulative distribution function (ECDF) of the sample data and the cumulative distribution function (CDF) of the reference distribution. It is defined as

[TABLE]

where sup is the supremum, $F$ and $\hat{F}$ denote the ECDF and CDF computed from the estimated probability density function. The D values and the $p$-values obtained using different models are given Table 11.

From Figure 6, we observe that the density of the normal distribution is too low near the middle, high in the midrange and quite low on the tails. While on the other hand, it is interesting to observe how all the stable models specifically the hybrid model approximate the data well over almost the whole range. Also from Table 11 the K-S goodness of fit test reveals that D is the smallest for the hybrid.

Tata Consultancy Services Limited (TCS.NS) and National Thermal Power Corporation Limited (NTPC.NS)

The daily returns of 1475 adjusted closing prices for the two components of Nifty, namely NTPC.NS and TCS.NS is obtained from Yahoo Finance for the time period, January 3, 2011 to December 31, 2016. Figure 7 shows the scatter and the contour plots of their returns respectively and reveals that the data is heavily skewed downwards and several points are away from the origin. Also from the contour plot it is clear that the distribution of the data is neither normal nor elliptical stable.

Thus, multivariate stable distributions are more appropriate than multivariate normal or elliptical distributions. The theoretical assessment of whether the data is multivariate normal is done using the MVN package in R and the test used is Henze-Zirkler test. The test reveals that the data is not multivariate normal. The details of this test on the bivariate data are shown in Table 12.

We now model the pair $\mathbf{X}=$( TCS.NS, NTPC.NS) using a bivariate stable distribution with $L=12$ points of masses for spectral measure given by \mathbf{t}_{l}=\mathbf{s}_{l}=\big{(}\cos(2\pi(l-1)/L,\sin(2\pi(l-1)/L\big{)}\in S^{2}, $l=1,2,\cdots,L$. The estimates of the location vector $\hat{\boldsymbol{\delta}}$ obtained using our proposed hybrid method are $\hat{\boldsymbol{\delta}}=(0.0003077732,-0.0001562171)$ while $\hat{\alpha}=1.85591$, obtained after taking the mean of $\hat{\alpha}_{1}=1.841935$ and $\hat{\alpha}_{2}=1.87007$. To fit a strictly stable distribution transform $\mathbf{X}$ to $\mathbf{X}-\hat{\boldsymbol{\delta}}$. The results of the estimates of the spectral measure after fitting a strictly stable distribution to the shifted data is shown in Table 13 and is compared with other well-known methods such as Mohammadi et al. (M) [22] and Teimouri et al. (T) [38]. We observe that the values of the estimated masses obtained via our method are very much closer to the other two methods.

6 Concluding Remarks

To conclude, we make the following observations in relation to our proposed method.

Our proposed hybrid method (univariate) is non-iterative in nature as no further improvement in the MSE’s of the estimators resulted after the first iteration. Also, as in the case of Kogon-Williams method, the regressions are performed only once using ordinary least squares as opposed to the numerous iterations for Koutrouvelis’ method where regressions are performed using generalized least squares. Thus, in terms of computational efficiency, our method is as good as the method of Kogon-Williams.

• 2.

Koutrouvelis makes use of look-up tables in order to perform regressions while Kogon-Williams eliminated the need to use look-up tables for estimation. Our proposed hybrid method shows that if we retain the look-up table and modify it in order to obtain the regression estimates using ordinary least squares, the performance of the estimates improves significantly in comparison to the estimates obtained using Kogon-Williams method in terms of accuracy and low MSE.

• 3.

The method of moments by Press is said to yield poor estimates, however, in our proposed hybrid method we have used this method to get the estimate of $\delta$ by suitably choosing two points at which the sample characteristic function is evaluated. Thus, by using the Press’ method we obtained a very good estimate of $\delta$.

• 4.

For the multivariate case, we make use of our proposed hybrid method to obtain the estimators $\hat{\alpha}$ and the shift vector $\hat{\boldsymbol{\delta}}$. The discretized estimator $\hat{\boldsymbol{\gamma}}$ of the spectral measure $\Gamma$ is obtained via the empirical characteristic method suggested by Nolan et al. [29] with slight modifications. Thus, in terms of computational efficiency and accuracy, the new method outperforms the method of Mohammadi et al. (M) [22] and Teimouri et al. (T) [38].

• 5.

Finally, we give two applications of our proposed method using financial data, where, the distribution of the datasets considered is stable. For the univariate data, K-S goodness of fit test shows that our method best fits the data in comparison to the other methods. Though the maximum likelihood method of estimation is said to give the most accurate estimate, however, it is computationally, the slowest when applied to the two financial data. For the bivariate data, we observe that the values of the estimated masses obtained via our method are very close to the values obtained through the method of Mohammadi et al. (M) [22] and Teimouri et al. (T) [38].