In this article, we introduce five statistical procedures in the field of Signal Analysis. They feature as practitioner’s exercises, mere snippets in fact, addressed at supplementing renowned statistical tests and select procedures pertaining to the received literature. In such guise, the procedures range from simple explanatory devices to novel algorithmic routines. The first procedure introduces and develops regime-switch testing applied to linear and nonlinear signals, whereas the second procedure provides a novel approach to efficient and automatic time-series smoothing. The third and the fourth procedures consider Bayesian analysis for the purpose of back/forecasting a vector dataset through a given time domain. The fifth procedure represents a novel method for computing the Augmented Dickey-Fuller stationarity test in a context of instability of both trend and subsample variances, which frequently characterize nonlinear signals.

Keywords: bayesian models, regime switching, stationarity, smoothing techniques, JEL classes: C1, C4

Signal Analysis currently features tests and procedures addressed at verifying certain hypotheses, such as linearity, stationarity, changes of regime, and the kind of information utilized in back/forecasting time series, possibly by calibration. In such context, methodological updates are required, especially when it comes to dealing with series characterized by high irregularities or when unobservable are involved as in the case of Kalman filtering.^1,2 Many Regime-Switching techniques, for instance, are affected by strict requirements such as stationarity and single-point shift^3,4 which require amelioration or updating. Similarly, time-varying parameters utilized for back/forecasting purposes are better computed in a Bayesian context, in order to exploit prior pieces of evidence based on their known distribution. Finally, current stationarity tests^5,6 in a nonlinear setting with dramatic changes in trend and subsample variances may require substantial corrections. Section 1 is devoted to defining random signals, whereas Section 2 tackles the issue of Regime-Switch testing. Section 3 produces a novel technique in the field of signal smoothing. The following two sections illustrate the steps for implementing Bayesian Monte Carlo simulation and calibration, whereas the final section is addressed at supplying corrected stationarity tests. Section 7 concludes.

Definition of signal

All of the aforementioned procedures share the assumption of a signal represented by a real-valued random

sequence which, for the discrete time notation $t\in \left[1,T\right]$ , $T<\infty ,$ may be expressed as follows

${\left\{X_t\right\}}_{t=1}^T\sim \text{ I}\text{.I}\text{.D}\text{.}\left(\overline{X},{\sigma }_X^2\right)$        (1)

With unknown distribution, mean $\overline{X}=T^{-1}{\displaystyle \sum _{t=1}^T{X}_t}$  and finite variance ${{\displaystyle \sigma }}_x^2<\infty$ . More specifically, Equation (1) may be decomposed Cramér⁷ in the following fashion

$X_t=X_t^{*}+S_t+{\epsilon }_t$ (2)

Where $X_t^{*}=\text{E}\left(X_t|{\Omega }_{t-1}\right)$  and ${\Omega }_{t-1}={\left\{X_{t-j}\right\}}_j^J$ is the information set available at $t-j$ ,$j\in \left[1,J\right]$ , $J\le T$ is a maximum preselect lag, and ${\epsilon }_t~\text{I}\text{.I}\text{.D}\text{.}\left(0,{\sigma }_{\epsilon }^2\right)$ . The first term $X_t^{*}$  of Equation (2) is the slow-mode component of the signal, either a linear trend or a smoother, or also a combination of the two. The second component $S_t$  is the periodical or quasi-periodical seasonal cycle, and the last component is a random Gaussian fast mode. Obviously, independent of $S_t$ , if the linear and/or the smoothed trends are not time-dependent, the first component zeroes out in mean and the process described by Equation (1) collapses to ${\epsilon }_t$. All components of $X_t$ are unobservable and must be retrieved from the raw dataset by means of appropriate computational methods. In particular the slow-mode term, if not solely represented by a linear trend, requires applying an efficient (i.e. minimum-variance) smoothing procedure. In addition, if $S_t$ is a high- or low-frequency phenomenon, it may be entrenched into either ${\epsilon }_t$ or $X_t^{*}$ and thus difficult to disentangle from either component. The signal so characterized constitutes a Monte Carlo Markov Chain⁸ where each observation of the signal $X_t$ may be unconditional or conditional upon a known distribution, abridged from prior information, or upon a transition probability ${\pi }_t\in \left[0,1\right]$.^9,10

Sequential and cluster-based regime switching

Regime switching is intended as a single or multiple change of pattern of the signal in terms of mean, variance, and autocorrelation structure.^9,10 It is often associated to rare events such as disasters with permanent effects like earthquake aftermaths or financial crises.¹¹ Changes associated to regime switching are generally abrupt and of sizable magnitude, and may be tested for their statistical significance by means of some preselect critical-value standard.¹² Regime switching is closely associated to structural-break testing.¹³ However, it is more pervasive by attempting to establish not only the date and size of breaks, but also the moment effects over the signal of interest. Moreover, regime switching is deeply involved with threshold behavior, which analyzes the event of a variable taking an extreme value at some time of its evolution. In practice, the parameterization of the signal after each break is tested for in a context of a Markov switching model. Many tests exist for the detection of regime switches,² whereby the different states of the signal may be identified within a given degree of likelihood. Among these tests, two stand out for their usefulness and simplicity:

1. Sequential F-test (two-sample Kolmogorov-Smirnov test) for linear (nonlinear) signals.^3,4
2. Sample cluster-based Quasi-Likelihood Ratio (QLR) test of the null hypothesis of one versus two regime switches Cho et al.,¹² or of m versus $m+1$ regime switches.¹⁴

The first kind of test in case of linearity (nonlinearity) is modeled on a local-to-global basis by comparing the variance (the tied-rank score sum) of subsamples to the variance (the mean tied-rank score sum) of the entire sample. Let ${\left\{x_{m,\tau }\right\}}_{\tau =1}^{{\rm T}}$ be the m.th subsample of Equation (1), for $m\in \left[1,M\right],\text{ 1}\le M\ll T$ and M the maximum number of subsamples, and also$\tau \in \left[1,{\rm T}\right],{\rm T}<T$. The subsample length is established arbitrarily, usually on a peak-to-peak basis. The F-test version for establishing the number of regime switches is given by the following statistic

$F_m=\frac{Var\left(x_{m,\tau }\right)}{Var\left(X_t\right)},\forall m\in M$ (3)

Where $Var\left(.\right)$ represents the variance of the braced argument. From Equation (3), the M-sized string of F statistics and associated p-values may be obtained. A regime switch, namely a subsample variance outlier with respect to the sample variance, is said to significantly exist if its m.th p-value falls short of some preselect critical value, usually 5%.

In case of a nonlinear signal, the test statistic corresponding to $F_m$ is represented by the two-sample Ansari-Bradley test for equal sample tied-rank scores, also defined as quasi-variances.¹⁵ The Ansari-Bradley test requires several steps departing from the padded signal vector $Z_{\theta }=\{x_{m,\tau };X_t\};\theta \text{=}\tau \text{+}t\in \Theta ;\Theta \text{=}{\rm T}\text{+}T\text{; }\forall m\in M$ , and from the associated vector $R_{\theta }$ of the locations of the values in $Z_{\theta }$ ranked in increasing order. Next, the locations of vector $R_{\theta }$ are parted into two group vectors depending on their origin, namely, their belonging to the subsample or to the entire sample. Formally, the concept may be expressed as $R_1=R_{\theta }\left|\theta =\tau \right.;R_2=R_{\theta }\left|\theta =t\right.$ . Subsequently, the tied-rank vector $U_{\theta }$ from each end of $Z_{\theta }$ is computed, so that the smallest and largest values of $U_{\theta }$ get rank 1, the next smallest and largest get rank 2, etc. Finally, the Ansari-Bradley test statistic¹⁵ is shown to be computed as follows

$W_m={\displaystyle \sum _{\tau =1}^{{\rm T}}{U}_{\tau }}\left|R_1\right.$ (4)

Which is the sum of the tied ranks associated to the m.th subsample? From Equation (4), the test statistic $W_m^{*}$ may be obtained by subtraction of the mean tied-rank score sum together with the corresponding p-value. The statistic is asymptotic normal under the null hypothesis that two stochastic processes come from the same distribution, against the alternative that they come from distributions sharing the same median but different variances. In the present context, it may also be loosely viewed as a test for heteroskedasticity.¹⁶ The QLR test, in its simplest form, checks the null hypothesis of one versus two regime switches, and is built up from forming two non-sequential clusters each characterized, possibly, by different distributions. Thereafter, the QLR is computed upon the two samples, with the proviso that these produce a Gaussian mixture where the standard Likelihood-Ratio testing cannot be implemented. For this purpose, the relevant literature provides appropriate critical values.¹² Because of its stricture in admitting a maximum of two regime switches and of its computational complexity, the QLR test receives here no applied treatment.

Spectral representation theorem for signal decomposition and optimal smoothing

A large amount of signal smoothers (i.e. smoothed trends) is nowadays available for applied work in many different scientific fields. Of these, the most notable are the Savitzky-Golay, the Daubechies, the Hodrick-Prescott filters and the Hilbert Huang Transform.^17–23 All of these methods, however, share serious problems when it comes to data resolution. In fact, they all require prior input of specific parameters, namely, the polynomial degree for the first two methods, the window width for the first and the decomposition level for the second method. In addition, it requires a preselected smoothing coefficient for the third method and, finally, a ‘sifting’ grid search for the last. For details of such inputs, the reader is redirected to the respective reference listed in this paper. Suffice here to warn that arbitrary, and thus casual, choice of any of these parameters is conducive to the risk of over/under fitting of the slow mode $X_t^{*}$ with respect to the original signal $X_t$ of Equation (1). The immediate consequence is inefficient estimation in terms of suboptimal denoising, namely and loosely stated, a bang-bang signal extraction yielding either a too flat smoother that is statistically inconsistent with the original peak/trough pattern or that too closely replicates the original signal itself. We propose here a fast, automatic, and efficient smoothing procedure that applies the Spectral Representation Theorem to a signal $X_1$ of whichever nature: stationary or not, linear or nonlinear. The signal may be approximated by means of the following periodic function

${\stackrel{⌢}{X}}_{t,k}=\mu +{\displaystyle \sum _{j=1}^J\left[{\varphi }_k\sin \left({\omega }_k\cdot \left(t-1\right)\right)+{\phi }_k\cos \left({\omega }_k\cdot \left(t-1\right)\right)\right]}$ (5)

where ${\widehat{X}}_{t,k}$ is a fitted sequence of the original signal, $\mu$ is the estimated mean of the fitted signal, $k\in \left[1,K\right]$ , $2\le K\ll T$ is a maximum lag integer, ${\left\{{\varphi }_k\right\}}_{k=1}^K\text{, }{\left\{{\phi }_k\right\}}_{k=1}^K$ are real-valued random coefficient sequences, both $\text{I}\text{.I}\text{.D}\text{.}$ with finite mean and variance. Finally, ${\left\{{\omega }_k\right\}}_{k=1}^K$ is a frequency sequence such that

${\omega }_k=\left\{\begin{array}{c}{T}^{-1}2\pi k\text{ if }\Delta X_t\text{=}f\left(X_{t-1}\right),\underset{k\to 0.5T}{Lim}\left({\omega }_k=\pi \right)\\ T^{-1}\pi k\text{ if }\Delta X_t\text{= }f\left(X_{t-1}^3\right)\text{, }\underset{k\to T}{Lim}\left({\omega }_k=\pi \right)\end{array}\right.$ (6)

Where $\Delta$ is a first-difference operator, $f\left(.\right)$ is a functional form, $\Delta X_t=f\left(.\right)$ expresses the existence of linearity and nonlinearity of the process as evidenced by the value of the exponent attached to the lagged signal. The fitted signal of Equation (5), which is the smoother estimable by Ordinary Least Squares (OLS), may be defined as the Spectral Representation Transform of the original signal with time-varying amplitude (Hamilton, 1994). If the signal is affected by more than one Regime Switches and we let $e_{t,k}=X_t-{\widehat{X}}_t$ where $e_{t,k}\sim \text{I}\text{.I}\text{.D}\text{.}\left(0,{\sigma }_{e,k}^2\right)$ , the Central Limit Theorem (CLT) is empirically found not to hold asymptotically as $k\to T$ . In fact, while for $k=1$ the smoother ${\hat{X}}_{t,k}$ converges towards the linear or nonlinear trend of $X_t$ , the CLT holds only for a limited range of $K\in \left[2,\alpha T\right)$ , for $\alpha \in \left(0.02,0.04\right)$ . In other words, for all lags $k\ge 2$ , $e_{t.k}^2$ tapers off, but then widens up again as K gets larger. Hence $L$ whereas $L$ for $1>\beta >\alpha$ , where the values of $\alpha$ are obtained by means of 1,000 Monte Carlo random normal and first-order autoregressive simulations of Equation (5). If the signal is not affected by Regime Switches, then the CLT applies and we have $\underset{k\to T}{\lim }\left({\stackrel{⌢}{X}}_{t,k}=X_t\right)$ . The optimal smoother is therefore defined as

${\widehat{X}}_{t,K^{*}}={\widehat{X}}_{t,k}\left|\min \left(P_k\right)\right.;\forall k\in \left(2\le K\right)<0.04T$ (7)

For the optimal lag $K^{*}$ of Equation (5) and the performance index closely akin to the Percent Root Mean Square Error (PRMSE), defined as

$P_k=\left(\frac{‖e_{t,k}‖}{‖X_t‖}\right)$       (8)

Where $||.||$ is the Euclidean norm of its argument.         

Bayesian monte carlo simulation with bootstrapping for back/forecasting

Bayesian Monte Carlo simulation with bootstrapping is one of the two signal back/forecasting procedures herein proposed. The other procedure is shown in Section 6. Whereas the latter is a full-blown Kalman filter with time-varying parameters extracted from Gibbs sampling,² this procedure is based on bootstrapping known dates and values of the signal to produce a feasible MCMC which represents the back/forecasted signal along select $P<\infty$ periods. This procedure is thus a game of chance which involves all of the elements of Equation (2) such that we may proceed to back/forecasting a signal by means of two different techniques:

1. Back/forecast the smoother of Equation (4) by P periods and thereafter append stationary shocks drawn from the distribution of the original signal,²⁴
2. Perform random reconstruction of the cyclical pattern of the original signal by applying the principles of peak/trough selection and bootstrapping these cyclical extrema without replacement.²⁵

While the first technique works well only under the assumption of fully or fairly regular cycles in amplitude and length, and is likely to be affected by substantial uncertainties,^26–28 the second technique preserves the stochastics of the original signal and is by consequence more trustworthy from the distributional viewpoint. There follows that Bayesian prior values drawn from the known distribution of the original signal can be utilized along the entire signal sample so as to forestall predictive uncertainties.

Implementation of this preferred technique requires three consecutive steps:

1. Find the peaks and troughs of the signal$X_t$  throughout the entire sample. Scrap all occurrences separated by a cyclical distance far exceeding or falling too short of the average cycle.
2. The number of trough-to-trough cycles remained may be denoted as $1\le C\ll T$ . Subsequently form a matrix $F\left(C,\text{3}\right)$  of triplets of dates in the following order: lower-date trough, peak date, and upper-date trough.

The latter is the lower-date trough of the subsequent cycle, and so on. Each triplet encapsulates a specific cycle of given length. Bootstrap N times the rows of matrix F so as to obtain a matrix $G\left(P,D\right)$  of the reshuffled values of the signal $X_t$ . For each row of this matrix such values might be further simulated up to D draws to produce Bayesian priors from the distribution of $X_t$ and generate a new matrix G with the posteriors. From this matrix find the series which mostly approximates the original signal $X_t$  in terms of mean and/or variance. This is the optimal prediction signal ${\widehat{X}}_p^{**}$ , for $p\in \left[1,P\right]$ , commonly defined as the “mode” since it is the closest to the actual data. This procedure is quick and easy to implement as opposed to the lengthier procedure featured in the next section, which needs a larger, often cumbersome informational set provided by the proxy instrument (s) utilized.

Bayesian calibration for back/forecasting

Calibration is a technique for back/forecasting the signal $X_t$ by means of a matrix of proxy instruments $Z_{\tau }:\left(\text{T,}n\right)$ , where $\tau \in \text{T>}T$ and $n\ge 1$ respectively are the sample length and size. It assumes a statistically significant relationship among the two datasets for the overlapping sample length T such that back/forecasting P periods ahead (behind) by means of $Z_p$ is expected to produce a mode ${\widehat{X}}_p^{**}$ . The relationship is customarily expressed in terms of a fixed-time parameter set $B_0:\left(n,1\right)$ usually estimated by means of OLS.²⁹ Bayesian calibration in the present context is an extension of the above model as it contemplates time-variable parameters connecting the two datasets and letting the prediction process be estimated by a dynamic state-space (SS) model whose upshot is a MCMC process characterized by ergodicity.² This model is represented by two equations whose workings in a Bayesian framework are better detailed elsewhere³⁰

$\begin{array}{l}{B}_{p+1}=B_p\left|{\Omega }_p+e_p\right.\\ X_p=C_p{B}_p+v_p\end{array}$ (9)

Where, for p given as above (Sect. 1.3), $B_p\text{ and }{X}_p$ respectively are the state(s) and the observable(s), whereas $e_p\text{ and }{v}_p$ are orthogonal I.I.D. zero-mean disturbances. More precisely, $B_p:\left(P,n\right)$ is the time-varying parameter set whose first input is given by the Bayesian prior $B_0$ conditional upon a normal distribution. The variable $C_p$ is a slope parameter set connecting the contemporary state(s) and observable(s) and has the same size as $B_p$ .

The real-time information matrix ${\Omega }_p=\left\{{\Pi }_p,Q_p,R_p\right\}$ includes the covariances ${\Pi }_p=B_p{B}_p\text{'},Q_p=e_p\text{'}{e}_p,R_p=v_p\text{'}{v}_p$ , all of which positive definite and of size $\left(n,n\right)$ except the last which is a scalar. The first component of ${\Omega }_p$ is the Riccati matrix of the state covariance, which is distributed as an Inverse Wishart, whereas the other two components respectively are the state and the observable error covariances, distributed as Inverse Gamma. Matrix ${\Omega }_p$ may be defined as the prior information set for predicting $B_{p+1}$ and its estimation errors $e_p$ . Once this value is computed conditional on the covariances and on their distributional properties, the ensuing posterior ${\Omega }_{p+1}$ is utilized for computing $B_{p+2}$ in the same fashion, and so on along the process. This prior-posterior-prior sequence of the information set is at the heart of the conditional adaptive Kalman Filter for which the conditional parameters of Equation (9) at any stage are optimally estimated and ends as soon as P is hit.1 Needless to say, optimality refers in the present context to the valuable information set ${\Omega }_p$ delivered stage-wise onto the parameters. More explicitly, let the Kalman Gain be

$H_p={\left(C_p{\Pi }_p{C}_p^{\text{'}}+R_p\right)}^{-1}$ (10)

Such that the Riccati matrix be represented overtime as the following nonlinear and recursive process

${\Pi }_p={\Pi }_{p-1}+Q_p-{\Pi }_p^{\text{'}}{C}_p^{\text{'}}{H}_p{C}_p{\Pi }_p$ (11)

Where by the sequential MCMC parameter series may obtain as follows

$B_{p+1}=B_p+{\Pi }_p{C}_p{H}_p{\nu }_p^{\text{'}}$ (12)

Needless to say, all of the elements in Equations (10) to (12) enter the real-time information matrix ${\Omega }_p$ . Optionally, in order to gain further information by letting $B_p\text{ and }{B}_{p+1}$ be conditioned directly by the observable $X_p$ , the covariance matrix $Q_p$ in Equation (11) may be replaced by the variance of the observable itself.

Proposed method for correction of stationarity test statistics

In this last Section a novel method is proposed for correcting the conventional stationarity test statistic in a setting of non-flat linear trends, and non-constant smoothers and cycles. Among such tests, the most notorious applicable to linear signals is the t-type Augmented Dickey-Fuller (ADF) test^5,6 which embeds an option for the intercept and/or for linear trend. Similarly, for nonlinear signals the Kapeitanos-Shin-Snell (KSS) t-type test is commonly applied with similar options.^31,32 The critical asymptotic values with preselected p-values of the test statistics are tabulated in the respective references. The null hypothesis implied by both test statistics involves first differencing of the signal, which in Econometrics is dealt with Autoregressive Conditional Heteroskedasticity (ARCH) testing in case of subsample changing variances.^33,34 We endeavor in the present context to follow a richer avenue, given the plethora of ad hoc tests and statistics available for deeply scrutinizing the performance of the signal in terms of changing trends, means and/or variances, and cyclical shapes. Among these we test, in a two-sample setting, for major changes that might have occurred within. Specifically, we test from Equation (2) the evolution of the slow-mode term $X_t^{*}$ , which includes a linear trend ${\tilde{X}}_t$ estimable by OLS, and the smoother ${\hat{X}}_{t,K^{*}}$ of Equation (7). We also test for equal variances by means of the F- or Ansari-Bradley test of Equations (3) and (4) respectively. We finally add a test for different subsample shape in terms of cyclical amplitude (e.g. Schwabe cycle).

Let the following: ${\Psi }_1={\tilde{X}}_T/{\tilde{X}}_1$ , ${\Psi }_2={\hat{X}}_{T,K^{*}}/{\hat{X}}_{1,K^{*}}$ , ${\Psi }_3={FFT{\left\{X_i\right\}}_{i=1}^N/FFT\left\{X_j\right\}}_{j=N+1}^T$ , where $FFT$ stands for the Fast-Fourier Transform of the associated signal and N is the length of the first subsample. Let also ${\Psi }_4$ be the p-value of $F_k$ or of the $W_k^{*}$ test statistic for $k=1$ , as from Equation (3) or (4). Under the null of constant trends, variances, and cyclical amplitudes all over the entire sample we have $\text{E}\left({\Psi }_j\right)=1;\forall j\in \left[1,4\right]$ , where $\text{E}\left(.\right)$ the expectation operator of the braced argument is. Under the alternative we have $\infty >{\Psi }_j>0;\forall j\in \left[1,4\right]$ . Under either or all the nulls, we have, for $\text{TS}=\text{ADF or KSS}$ of the raw dataset, $\text{E}\left(\text{TS}-{\displaystyle \sum _{j=1}^4\text{E}\left({\Psi }_j\right)}-4\right)=\text{TS}$ , such that the corrected stationarity test statistic is

${\text{TS}}^{*}=-\left(\left|\text{TS}\right|-\left|{\displaystyle \sum _{j=1}^4\left|{\Psi }_j-1\right|}\right|\right)$ (13)

Where the computed TS are penalized whenever one or more of the coefficients ${\Psi }_j$ exceeds or falls short of unity. Penalization consists of the reduction in the absolute value of the computed TS, possibly below the associated critical value or even turning into a positive value. In practice, the farther away are the coefficients ${\Psi }_j$ from their expected values, the greater the risk of nonstationarity in spite of stationarity originally detected by the conventional test statistics. Incidentally, the asymptotic 10%, 5%, and 1% critical values of the computed TS are -2.22, -2.93, and -3.40.^35–37 Incidentally, such values are close to some simulated experiments over 10,000 replications of the nonlinear version of Equation (2).

Five statistical practitioner’s procedures in support of signal analysis have been introduced and discussed. Such snippets were shown to be useful for better and efficient implementation of the corresponding procedures and of some select hypothesis tests existing in the received literature. Of notable consideration are the optimal smoothing and the corrected ADF test procedures, both characterized by significant novelty in their implementation techniques. The Bayesian stepwise procedures for Monte Carlo simulation and calibration are also of interest because of their ability to perform back/forecasting of signals over long periods of time, if necessary. Of the two, the former may be preferred because quick and easy to implement. The first and the last snippets were found to be particularly useful in a highly erratic environment regarding the signal. The latter, in particular, is worth of notice because of its capability to reverse the results of the traditional stationarity test statistics.

None.

1. Kalman RE. A new approach to linear filtering and prediction problems, Trans. ASME. 1960;82:35–45.
2. Carter CK, Cohn P. On Gibbs sampling for state space models, Biometrika. 1994;81:541–553.
3. Rodionov SN. A brief review of the regime shift detection methods, manuscript, Joint Institute for the Study of the Atmosphere and Ocean. University of Washington, Seattle, WA. 2005.
4. Rodionov SN, Overland JE. Application of a sequential regime shift detection method to the Bering Sea ecosystem. ICES Journal of Marine Sciences. 2005;62(3):328–332.
5. Dickey DA, Fuller WA. Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association. 1979;74(36):427–431.
6. Said E, Dickey DA. Testing for unit roots in autoregressive moving average models of unknown order. Biometrika. 1984;71(3):599–607.
7. Cramér H. On some classes of nonstationary stochastic processes, Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability. Berkeley and Los Angeles, University of Los Angeles Press. 1961;57–77.
8. Asmussen S, Glynn PW. Stochastic simulation: algorithms and analysis, Springer Science and Business Media. 2017;ISBN: 978-0-387-30679-7, pp. 470–476.
9. Hamilton JD. Time Series Analysis, Princeton University Press, Princeton, N.J., 1994. 1994. p. 470–799.
10. Hamilton JD. Regime-switching models, New Palgrave Dictionary of Economics, 2nd edition, Durlauf S, et al., editors. Palgrave McMillan Ltd., 2008. 2007.
11. Ang A, Timmermann A. Regime changes and financial markets. Annual Review of Financial Economics Annual Reviews. 2012;4:313–337.
12. Cho JS, White H. Testing for Regime Switching. Econometrica. 2007;75(6):1671–1720.
13. Bai J, Perron P. Estimating and testing linear models with multiple structural changes. Econometrica.1998;66(1):47–78.
14. Kasahara H, Shimotsu K. Testing the Number of Regimes in Markov Regime Switching Models. Mimeo London School of Economics. UK. 2017. p. 1–68.
15. Ansari AR, Bradley RA. Rank-sum tests for dispersion. Annals of Mathematical Statistics. 1960;31(4):1174–1189.
16. Breusch TS, Pagan AR. A simple test for heteroscedasticity and random coefficient variation. Econometrica. 1979;47(5):1287–1294. 
17. Savitzky A, Golay MJE. Smoothing and differentiation of data by simplified least squares procedures, Analytical Chemistry. 1964;36(8):1627–1639.
18. Daubechies I. Ten lectures on wavelets, SIAM, Society for Industrial and Applied Mathematics. 1992. p.1–357.
19. Hodrick R, Prescott EC. Postwar U.S. business cycles: an empirical investigation. Journal of Money, Credit, and Banking. 1997;29: 1–16.
20. Ravn MO, Uhlig H. On adjusting the Hodrick-Prescott filter for the frequency of observations. Review of Economics and Statistics. 2002;84(2):371–376.
21. Huang NE, Shen Z, Long SR, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society in London. 1998;454:903–993.
22. Huang NE, Wu ML, Qu W, et al. Applications of Hilbert–Huang transform to non-stationary financial time series analysis. Applied Stochastic Models in Business and Industry. 2003;19(3):245–268.
23. Huang NE, Wu Z. A review on Hilbert-Huang transform: Method and its applications to geophysical studies. Review of Geophysics. 2008;46:RG2006.
24. Gamerman D, Lopes HF. Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference, 2nd edition, Chapman & Hall/CRC Press. 2006. pp. 344.
25. Travaglini G. An econometric investigation of the sunspot number record since the year 1700 and its prediction into the 22nd century. Advances in Space Research. 2015;56(5):992–1002.
26. Mann ME, Bradley RS, Hughes MK. Northern Hemisphere temperature during the past millennium: inferences, uncertainties and limitations. Geophysical Research Letters. 1999;26(6):759–762.
27. Mann ME, Jones PD. Global surface temperatures over the past two millennia. Geophysical Research Letters. 2003;30(1):1880–1881.
28. Von Storch H, Zorita E, Jones JM, et al. Reconstructing past climate from noisy data. Science. 2004;306(5696):679–682.
29. Mann ME, Zhang Z, Hughes MK, et al. Proxy-based reconstructions of hemispheric and global surface temperature variations over the past two millennia. Proc Natl Acad Sci U S A. 2008;105(36):13252–13257. 
30. Travaglini G. Testing the hockey-stick hypothesis by statistical analyses of a large dataset of proxy records. Pattern Recognition in Physics. 2014;2(2):35–57.
31. Kapeitanos G, Shin Y, Snell A. Testing for a unit root in the non-linear STAR framework. Journal of Econometrics. 2003;112:359–379.
32. Kapeitanos G, Shin Y. Testing the null hypothesis of nonstationary long memory against the alternative hypothesis of a nonlinear ergodic model. Econometrics Review. 2011;30(6):620–645.
33. Engle RF. Autoregressive Conditional Heteroscedasticity with estimates of the variance of the United Kingdom inflation. Econometrica. 1982;50:987–1007.
34. Bollerslev T. Generalized Autoregressive Conditional Heteroskedasticity. Journal of Econometrics. 1986;31(3):307–327.
35. Harvey DI, Leybourne SJ, Xiao B. A powerful test for linearity when the order of integration is unknown. Studies in Nonlinear Dynamics & Econometrics. 2008;12(3):1–2.
36. Kolmogorov AN. Sulla determinazione empirica di una legge di distribuzione. Giornale dell’ Istituto Italiano degli Attuari. 1993;483–491.
37. Massey FJ. The Kolmogorov-Smirnov test for goodness of fit. Journal of the American Statistical Association. 1951;46(253):68–78.