Random coefficient autoregressive processes and the PUCK model with fluctuating potential

The potentials of unbalanced complex kinetics (PUCK) model consists of a random walker subjected to a potential centered at its moving average position. We study the PUCK model with fluctuating quadratic potential, showing that it is a special case of the random coefficient autoregressive (RCAR) process and thus a member of the same class of processes as the autoregressive conditional heteroskedasticity (ARCH) process; both RCAR processes but with different coefficient dependence. For the general Gaussian RCAR process, we use generalized Fibonacci numbers to derive explicit expressions for the mean squared displacement and the autocovariance function. We also obtain the conditions for divergence of variance, which imply tail index of the power-law-tailed cumulative probability distribution. Translating the results to the PUCK model, we analyze US Dollar/Japanese Yen market data and demonstrate that the model simultaneously reproduces empirical facts of price time series: diffusion and autocorrelation of prices and heavy-tailed distribution of price changes.


Introduction
Stochastic processes are utilized as mathematical models for the time evolution of systems presenting hard-to-track complicated dynamics, such as the Brownian motion of particles suspended in a fluid, whose positions at each instant of time are then represented by random variables [1,2]. One of the basic stochastic processes used to model time series is the autoregressive (AR) process, defined as: i.e. the present value of the process is determined by a linear combination of the past n values, with coecients φ u real numbers, and an error term ξ(t) which is part of a sequence of independent and identically distributed random variables. Some of the first applications of the AR process in the description of real time series can be credited to Yule and Walker with the analysis of sunspot numbers and atmospheric pressure variations, respectively [3,4]. Since then many refinements and generalizations have been proposed, such as the composition with the moving-average https://doi.org/10.1088/1742-5468/aaf109 J. Stat. Mech. (2019) 013403 process-linear combination of error terms-from which originated the autoregressivemoving-average (ARMA) process and the integration of the ARMA process leading to the autoregressive-integrated-moving-average (ARIMA) process [5]. Those models have become a standard basic framework for the characterization and forecasting of time series, being used for the analysis of real data in diverse fields, e.g. biology and medical science [6][7][8], climate science [9][10][11], hydrology [12,13], study of trac [14], and tourism management [15].
A particular area in which the AR process has marked importance is in econometrics for the modeling of financial time series [16][17][18][19]. However, the simple linear AR process fails in capturing some empirical facts of market price time series, in special the heavy-tailed probability distribution of price changes and volatility clustering (according to Mandelbrot: 'Large changes tend to be followed by large changes-of either sign-and small changes tend to be followed by small changes') [20][21][22]. An extension of the AR process proposed to account for these observations in financial time series was introduced by Engle [23,24], called the autoregressive conditional heteroskedasticity (ARCH) process: In the ARCH process, it is the conditional variance of the present value of the process that is determined by a linear function of the past n values. In fact, the ARCH process can be viewed as a member of the wider class of random coecient autoregressive (RCAR) processes, proposed some years before ARCH as a development of the linear AR process outside the context of econometrics [25,26]: where η j (t) are taken here as identically distributed but not necessarily independent. Generally, the ARCH process and the corresponding RCAR process, i.e. RCAR with null coecients φ j and mutually uncorrelated η j (t), are second-order equivalent, denoting that they have the same conditional mean and conditional variance [27] (note that the ARCH process defined in equation (2) corresponds to the case ARCH(0,q) of the generalized form ARCH(p,q) considered in [27]); the identification between the two processes is exact when the random variables are all normally distributed. Besides the varying conditional variance, a remarkable property of RCAR processes is precisely the power-law-tailed distribution even when η j (t) are light-tailed-distributed [28][29][30].
Still in the line of models for financial time series, now in the recent field of econophysics, another process developed to represent market price dynamics is the potentials of unbalanced complex kinetics (PUCK) model, which describes a random walker subjected to a potential with its center at the moving average of its last positions [31][32][33][34]: where U n (x, t) is the possibly time-dependent potential and x n (t) = 1 mainly focused on the time-constant quadratic potential, although higher-order potentials were also considered to deal with financial crisis periods or interventions in the market [35][36][37]. The time-constant quadratic potential corresponds to a linear force, with the derived PUCK model having the same structure as a linear AR process and as such not producing a heavy-tailed distribution.
In this work, we extend the time-constant quadratic potential by adding random fluctuations to it, resulting in a (time-varying) fluctuating quadratic potential. We show that the PUCK model with fluctuating quadratic potential is also an RCAR process having a similar structure to the ARCH process but with distinct dependence relations among the coecients. Such identification with the RCAR process motivates us to first study the general RCAR process and then translate the obtained theoretical results to the PUCK model. We concentrate on three properties observed in real financial time series: deviations of the mean square displacement from the ordinary diusion for short time scales, the quick decay of the autocorrelation function, and the already mentioned heavy-tailed distribution [31,[38][39][40][41]. By analyzing real data from the foreign exchange market, we obtain quantitative agreement between the PUCK model with fluctuating quadratic potential and those empirical facts of the market.
In the next section, we present the PUCK model with fluctuating quadratic potential and its RCAR structure, motivating the study of general RCAR processes. We start the investigation with the ordinary AR process in section 3, deriving the explicit expressions of the mean squared displacement (MSD) and the autocovariance function (and consequently the autocorrelation function). We then use those results in section 4, where we study the Gaussian RCAR process, finding that the MSD and the autocovariance have the same form as the AR process, just amplified by the random fluctuations. We also discuss the power-law tail of the cumulative probability distribution and the conditions for tail index α = 2. Finally, we translate the theoretical results of the general RCAR process to the PUCK model in section 5 and, in section 6, we analyze real market data using the PUCK model with fluctuating quadratic potential and show the reproduction of the market empirical facts.

PUCK model with fluctuating quadratic potential
The PUCK model defined in equation (4) is characterized by the potential U n (x, t). From Taylor expansion, the first approximation for an arbitrary time-constant potential is the quadratic potential [31][32][33][34]: (5) Figure 1 shows simulations of the PUCK model for n = 5 and Gaussian noise ξ ∼ N (0, 1) together with the representation of the associated quadratic potentials. The case b = 0 (black line) corresponds to the absence of potential (or to a constant potential), producing an ordinary random walk with all steps independent. Case b > 0 (violet and blue lines) refers to a stable potential: the walker is attracted to the average of its previous positions and the process is less diusive than a random walk. Case https://doi.org/10.1088/1742-5468/aaf109 J. Stat. Mech. (2019) 013403 b < 0 (orange and red lines) refers to an unstable potential, repelling the walker from the average of its previous positions and making the process diusion faster than a random walk.
The concepts of stable and unstable potentials are particularly appealing for the price movement description, with the potential reflecting the conditions of the market and indicating the risks of transaction. From a microscopic point of view, it has been suggested that the potentials are constructed by the compositions of psychology of the traders and its stability is determined by their dominant strategy when facing market price changes: a stable potential arises when traders are contrarians (opposed to trends) and an unstable potential appears when they are trend-followers [42,43].
Although able to reproduce diusion and correlation properties of price dynamics, the PUCK model with time-constant quadratic potential (equivalent to a linear AR process) does not produce a price changes heavy-tailed distribution for light-tailed distributed noise. Also, considering that the potential is built from the psychology of traders, the existence of some fluctuations related to their degree of strategy flexibility is expected. A previous work added stochasticity to the potential function by assuming that the parameter b follows an Ornstein-Ulenbeck process [34], which can be seen as an RCAR process with coecients correlated in time [44]. Such time-evolving potential yields heavy-tailed distributions but it restricts fluctuations around zero and makes the derivation of analytical results dicult. Following the idea of adding random fluctuations to the potential but aiming at an analytically tractable process with straightforward interpretation, here we introduce the fluctuating quadratic potential: where θη(t) is the random fluctuation in the potential associated with the flexibility of traders in the market (or thermal fluctuation if we are considering physical potentials). Of course, θ = 0 recovers the time-constant potential case. Inserting the potential (6) in the definition of the PUCK model (equation (4)), we verify that the PUCK model with fluctuating quadratic potential is a special case of RCAR process (compare with equation (3)): Dierent from the ARCH process, which is defined from the conditional variance, the exact correspondence between the PUCK model with fluctuating quadratic potential and the RCAR process with η j (t) = η(t), ∀j, is direct, without requiring Gaussian random variables. Observe that, even though the PUCK model with fluctuating quadratic potential and the ARCH process are members of the same class of RCAR processes, there are no specifications of this PUCK model equivalent to an ARCH process because of the distinct coecient dependence in each case. To study the properties of the PUCK model with fluctuating quadratic potential, in the next sections we consider the general RCAR process and investigate its diusion and correlation properties and its probability distribution.

AR process
We start with the ordinary AR process, which is the limit case of the RCAR process without randomness in the coecients. Results obtained in this section are used in the next one on general RCAR processes.
The AR process is given by the linear nonhomogeneous dierence equation of order n represented in equation (1). In dierence equations theory [45] and time series analysis [5,[16][17][18][19] one usually employs the backshift operator representation Bx(t) = x(t − 1). We choose an alternative approach, using the generalized Fibonacci numbers F (t) [46][47][48]: note that the recursive expression for F (t) represents the corresponding homogeneous dierence equation with Fibonacci initial conditions: F (t < n − 1) = 0, F (n − 1) = 1. This approach is useful when generalizing to RCAR processes (section 3). By recursion, we express the causal solution of the AR process with the generalized Fibonacci numbers as: The graph on the right depicts the potential functions U 5 (x): b = 0 represents the absence of potential yielding an ordinary random walk; b > 0 (b = 5 and b = 2) gives a stable potential resulting in processes less diusive than the random walk; and b < 0 (b = −1 and b = −1.5) refers to an unstable potential with a corresponding process more diusive than the random walk. https://doi.org/10.1088/1742-5468/aaf109 An expression for F (t) can be obtained (appendix A) and it is given in terms of the r n distinct roots λ j , each with multiplicity m j , of the characteristic polynomial p(λ) = λ n − n j=1 φ j λ n−j : with: Note that an AR process of order n cannot have a zero-root λ j = 0, otherwise it would be of order n − 1.
The AR solution representation in equation (9) and the expression of F (t) in equation (10) allow us to compute the MSD . We consider the case of a single-unit root λ 1 = 1, m 1 = 1, |λ j =1 | < 1 and take ξ(t) = ση 0 (t), η 0 ∼ N (0, 1), i.e. Gaussian noise with variance σ 2 . The expression of MSD of AR in this case is: where p * (λ) = 1 − n j=1 φ j λ j is the reciprocal polynomial of p(λ). For large time scale τ , the linear term due to the unit root λ 1 = 1 is dominant and the process has the same behavior as the ordinary diusion of a random walk; for short time scale τ , the For stationary AR processes, with |λ j | < 1, ∀j, the MSD does not contain the linear term and when τ → ∞ it converges to twice the variance of the process. The autocovariance function C(τ ) in this case is (see appendix B): The stationary AR process with Gaussian noise is a Gaussian process and therefore the position x is normally distributed with zero mean and variance C AR (0). In the following section, with the introduction of random coecients, this Gaussian distribution is replaced by a power-law-tailed distribution.

RCAR process
The RCAR process is defined by equation (3), which is classified as a linear nonhomogeneous dierence equation of order n with variable coecients, in this case, random coecients. Analogous to the previous section on AR processes, calling we define the generalized random Fibonacci numbers F s (t): A representation of the causal solution of the RCAR process using F s (t) is obtained by recursion (in a similar way to [49]): which resembles equation (9) for the AR process, but now coecients φ v (t) are timedependent and the index s of the generalized random Fibonacci numbers indicates the time shift of the coecients in the definition of F s (t).
An expression of the proposed generalized random Fibonacci numbers can be written in a convenient way using the deterministic generalized Fibonacci numbers naturally depending on the whole history of the time-dependent fluctuations η v (t).
Then we compute the MSD D 2 (τ ) of the RCAR process, considering again Gaussian noise ξ(t) = ση 0 (t), η v ∼ N (0, 1) and Gaussian fluctuations of the coecients η v ∼ N (0, 1), ∀v; η 0 (t) and each η v =0 (t) are independent but there can be dependences among η v =0 (t). Either for the case of a single-unit root or the stationary case, we have (appendix D): i.e. the MSD for an RCAR process has the same functional form as the one for the corre sponding AR process with the same coecients φ v but amplified by a factor 1 1−γ , provided that |γ| < 1. Similarly, the autocovariance function C(τ ) of a stationary RCAR process is: which recovers the result that the autocorrelation function C RCAR (τ ) C RCAR (0) of a stationary RCAR process is the same as the corresponding AR process [25].
The amplification factor 1 1−γ depends on the intensities θ v of the fluctuations of the coecients and γ is given by (appendix D): https://doi.org/10.1088/1742-5468/aaf109 which is always non-negative.
For an explicit expression of γ we consider two situations: it is only non-zero when t = n − 1 and then: Observe that the ARCH process (equation (2)) corresponds to this case when all fluctuations η v are independent and the result agrees with the known variance of the ARCH process: 2. Non-zero roots of the characteristic polynomial p(λ): λ j = 0, ∀j In this case we use the expression of the generalized Fibonacci numbers in equation (10) and obtain: where q(λ) = n j=1 θ j η j λ j is the random polynomial containing the information about the fluctuations. Dierent from the zero mean coecients case, here γ is influenced by the dependence among η v . We study two extreme cases corresponding to the dependences in the PUCK and ARCH processes and calculate the term q(λ)q 1 λ : Note that for RCAR processes with ARCH dependence γ always diverges if there is a unit root λ 1 = 1 because p * (1) = 0 and q(λ)q 1 λ > 0 (considering θ j = 0, ∀j). This divergence does not happen for PUCK dependence provided that n j=1 θ j = 0 and consequently q(λ 1 )q 1 It is known that stationary RCAR processes have probability distributions with a power-law tail [30], a result derived using Kesten's theorem on the general matrix stochastic recurrence equation X(t) = A(t)X(t − 1) + B(t), with A(t) a real random matrix and B(t) a real random vector [50][51][52]. The main result for the RCAR process under study is that the (complementary) cumulative distribution function P (X > |x|) has a power-law tail with the tail index determined by: where . denotes a matrix norm and the random matrix Φ(s) is: In the simplest case of an RCAR process of order n = 1, equation (24) reduces to: Aside from the one-dimensional case, the multiplication of random matrices in equation (24) makes the derivation of a general explicit expression dicult, so that currently the relationship between RCAR parameters and tail index α can only be estimated via simulations. Nevertheless, we can obtain here the explicit expression for the special case α = 2: since from Kesten's theorem it is known that the probability distribution has a power-law tail, the divergence of the second moment x 2 gives the conditions for the cumulative distribution function tail index α = 2; this corresponds to the divergence of the variance of the process, which first occurs when γ = 1. Thus, γ = 1 implies α = 2.

Theoretical results for PUCK model
After studying the general RCAR process, we can now apply the derived results to the PUCK model with fluctuating quadratic potential (equation (7)). We first observe that the characteristic polynomial of PUCK p(λ) = λ n − 1 − b n λ n−1 − b n(n−1) n j=2 λ n−j always has a unit root λ 1 = 1 and the other roots obey |λ j =1 | < 1 for −2 < b < 2(n − 1) if n is even or −2 < b < 2n if n is odd. Also, all roots are taken as having multiplicity m = 1 (numerically tested in [32]). Then, within the limits for parameter b and taking Gaussian random variables, we can use equations (12) and (17) valid for γ < 1 and with γ, noting that the sum of fluctuation intensities − θ n + n j=2 θ n(n−1) = 0: where q(λ) = − θ n λ + θ n(n−1) n j=2 λ j . The exact diusion behavior represents a new result for the PUCK model, even for the time-constant quadratic potential, with previous studies using approximations [34]. Figure 4 shows the square root of the MSD D(τ ) for the PUCK model of order n = 5. In figure 4(a), we consider the PUCK model with time-constant quadratic potential (θ = 0) and confirm the deviations from the ordinary diusion D(τ ) ∼ τ 1 2 at short time scales when b = 0. For b > 0 (violet and blue symbols/lines) the diusion is slower than a random walk at short time scales and for b < 0 (orange and red line) it is faster; for b = 0 (black symbols/line) the process is a random walk and presents an ordinary diusion. In figure 4(b), we examine the eect on the MSD of the fluctuation intensity θ in the PUCK model with fluctuating quadratic potential and parameter b = −1. The functional form of D(τ ) remains the same as the θ = 0 case (black symbols/line) but amplified by the fluctuations; we note that as θ approaches the divergence limit γ = 1 it becomes dicult to simulate extreme values and for the simulation to reach the theoretical value (blue line).
This PUCK model (of order n) is not stationary due to the presence of a unit root. Thus we consider its increments by dierencing the original process and obtain a stationary RCAR of order n − 1 (respecting the limits for parameter b and condition γ < 1): the characteristic polynomial for the increments p ∆ (λ) has the same roots as p(λ) with the exception of the unit root. From equations (13) and (18), the autocovariance function of the stationary process is then: where q ∆ (λ) = − θ n(n−1)

J. Stat. Mech. (2019) 013403
In figure 5, we depict the autocovariance function C(τ ) for the increments ∆x of the PUCK model using the same specifications as figure 4, considering time-constant quadratic potential (figure 5(a)) and fluctuating quadratic potential with b = −1 ( figure  5(b)). Negative values of the autocovariance occur for b > 0-stable potential attracts the walker to the average of its previous positions-while for b < 0, the autocovariance function remains non-negative. Equal to the MSD, the eect of fluctuations is only an amplification of the function C(τ ).
As discussed in the previous section for stationary RCAR processes, the increments of the PUCK model present a power-law-tailed distribution whose tail index is specified by equation (24). The value of parameter θ that produces tail index α = 2 can be exactly determined by analyzing the divergence condition of the variance C PUCK (0), which happens when γ ∆ = 1. Figure 6 illustrates the cumulative distribution function P (X > |∆x|) of the absolute increments |∆x| of the simulated PUCK model of order n = 5 with fluctuating quadratic potential and parameter b = −1. We cover several values of θ, from θ = 0 (black symbols) giving a Gaussian distribution, to θ = 2.706 which yields a power-law-tailed distribution with α = 2.
Finally, we construct in figure 7 the curves γ ∆ = 1 relating the parameter b and θ for a given order n that produces a power-law-tailed distribution with tail index α = 2

J. Stat. Mech. (2019) 013403
of the increments of the PUCK model with fluctuating potential. Although we do not have the exact curves representing other tail indices for n > 2, a curve γ ∆ = 1 divides the parameter space in two regions: below the curve are the parameter pairs that result in tail indices α > 2, and above the curve tail-indices α < 2.

PUCK analysis of market data
In this last part, we employ the PUCK model with fluctuating quadratic potential to analyze foreign exchange market data from the Electronic Broking Service (EBS). The utilized dataset contains the record of the traders' buy and sell quotes, i.e. desired transaction prices, with a 0.1 s resolution. We define the mid-quote as the average of the highest buy and the lowest sell quotes and construct its time series using tick time, i.e. a value is registered only when there is a change in the mid-quote. We analyze the mid-quote time series for the currency pair US dollar (USD) and Japanese yen (JPY) on August 30 2011 from 0:00 GMT (9:00 JST) to 8:00 GMT (17:00 JST); no major event occurred during this period but a highlight of August 2011 was the intervention by the Japanese government on August 4 contributing to a stabilization of the yen [53]. The mid-quote time series is presented in figure 8 together with the corresponding MSD, cumulative distribution function, and the autocorrelation function for changes in the mid-quote. We also show the autocorrelation function for squared mid-quote changes, which we have not studied analytically but it is also related to a stylized fact of market price time series: the long-range correlations of the magnitude (or squares) of price changes, a manifestation of the volatility clustering [38].
Parameter estimation procedures for the general RCAR exist, with least squares estimators and maximum likelihood estimators [54,55]. Here we use the commonly employed Yule-Walker equations to estimate the parameter b of the PUCK model from the autocorrelation function ρ PUCK (τ ) = C PUCK (τ ) C PUCK (0) for increments and exploit the results from the previous section to obtain the parameter θ directly from the tail index α of the cumulative distribution function and the parameter σ from the variance C PUCK (0). First, we consider the increments of the PUCK model (equation (29)) and select n = 2 from the functional form of the autocorrelation function for the mid-quote changes ( figure 8(e)), resulting in a stationary RCAR process of order 1. We also extract from the data the sample autocorrelation of lag 1 and the sample variance with corresponding standard errors: ρ PUCK (1) = −0.15 ± 0.01 and Ĉ PUCK (0) = (5.2 ± 0.1) × 10 −6 . From the Yule-Walker procedure, we have that − b 2 = ρ PUCK (1) and then the estimate b = 0.30 ± 0.02, where the uncertainty is linearly propagated from the standard error of the sample autocorrelation. Next, we can use the condition γ ∆ = 1 (that gives tail index α = 2) to find an upper bound for parameter θ. But since n = 2 for the analyzed market data, θ can be explicitly estimated through equation (26) for a Gaussian random variable:  (4 − b 2 − θ 2 )C PUCK (0) and the estimate is σ = 0.0016 ± 0.0001. In figure 9, we simulate the PUCK model with fluctuating quadratic potential using order n = 2 and parameters b = 0.3, θ = 1.4, and σ = 0.0016 and compare its statistical properties with the ones of the USD/JPY mid-quote time series. Having estimated the parameters directly from the statistical properties, we observe quantitative agreement for the MSD (figure 9(c)), the cumulative distribution function (figure 9(d)), and the autocorrelation function for increments (figure 9(e)). Regarding the autocorrelation function for squared increments, although dicult to determine if it follows a powerlaw for the studied dataset-as obtained when modeling using the autoregressive fractionally integrated moving average process (ARFIMA) [56,57]-it is clear that the PUCK model with fluctuating quadratic potential as we defined does not reproduce it ( figure 9(f)).
The disagreement in the autocorrelation function for squared increments suggests the existence of some kind of time correlation in the fluctuations of the potential (6). Based on the results in [34], with no particular empirical motivation but as a mathematical extension, we test the potential: where:   that is, the fluctuation β(t) follows an AR process of order 1, resulting in an AR process within another AR process. Since for this potential we do not yet have the exact expression of the relevant statistical properties, we perform numerical experiments and find that an autocorrelation function for squared increments similar to the one observed in the USD/JPY mid-quote time series is produced when the AR process defining the potential fluctuation (34) has a characteristic polynomial with root close to 1 (unit root). However, the introduction of time correlations in the fluctuation changes the functional form of the MSD and autocorrelation function for increments (with respect to equations (27) and (30), respectively) and the nature of the tail of the cumulative distribution function is obscured. Figure 10 depicts the PUCK model with the potential (33) using n = 2, b = 0.3, θ 1 = 0.999, θ = 0.0185, and σ = 0.002 15, and the comparison between its statistical properties and the ones of the market data. Detailed investigations on this type of process still need to be made to verify if it can reproduce (and how) those four characteristics of market price time series simultaneously and if it is actually justified from the perspective of behavioral economics.

Final remarks
Summarizing this work, we proposed the PUCK model with fluctuating quadratic potential as an extension of the PUCK with time-constant quadratic potential, the latter being able to reproduce diusion and some autocorrelation behaviors observed in market price time series but not producing a heavy-tailed probability distribution for the increments. The observation that the PUCK model with fluctuating quadratic potential is in fact a special case of RCAR process motivated us to conduct a theoretical study of the general RCAR, obtaining results that can be applied to other models such as the ARCH process.
We started with the simple AR process, deriving the explicit expression for the MSD and autocovariance function using generalized Fibonacci numbers, and then extended these results to the RCAR process, finding that the MSD and autocovariance function have the same functional form as the deterministic coecients case but with an amplification dictated by a factor γ. We recalled the known fact that RCAR processes have a cumulative distribution function with a power-law tail; the determination of the tail index α of this distribution is not trivial but we were able to find the condition for the particular case of α = 2: γ = 1.
Finally, we returned to the PUCK model with fluctuating quadratic potential, writing explicitly its MSD and the autocovariance function for increments and establishing the relations between parameters b and θ that yield a tail index α = 2. The appropriateness of the PUCK model with fluctuating quadratic potential was tested by analyzing foreign exchange data, the mid-quote time series of the currency pair USD/JPY on August 30 2011. Quantitative agreement was obtained for the MSD, the cumulative distribution function, and the autocorrelation function for increments, but the PUCK model with fluctuating quadratic potential (or any other general finite-order RCAR) was not able to replicate the empirical autocorrelation function for squared price changes, which is another important characteristic of price time series. The introduction of temporal correlations in the fluctuations of the quadratic potential, e.g. an AR structure, appears as a direction for future studies. https://doi.org/10.1088/1742-5468/aaf109

J. Stat. Mech. (2019) 013403
The PUCK model with fluctuating quadratic potential presents a reduced number of parameters (2+1 parameters) when compared to general RCAR (2n + 1 parameters) but it can still reproduce some of the statistical properties of market data, i.e. even though it represents a much smaller subset of possible RCAR processes, it has the necessary structural features for the modeling of the price time series. Some of those features are the presence of a single unit root while keeping n j=1 θ j = 0, allowing at the same time a linear MSD for large time scales and heavy-tailed distribution, and the specific dependence among the coecients that preserve the dierenced PUCK as a (stationary) RCAR. Other than that, the interpretation in terms of a potential reflecting the psychology of the traders is interesting to characterize the market condition and opens the possibility of considering other potential functions, although the analytical study of higher-order potentials is challenging due to the non-linear character of the corresponding processes.
where F(t) is a vector containing generalized Fibonacci numbers: 2) and Φ is the companion matrix of the polynomial p(λ) = λ n − n j=1 φ j λ n−j : By recursion, we can write: There are m j linearly independent generalized eigenvectors v jk , 1 k m j , corresponding to the eigenvalue λ j [58]. They satisfy: where 1 is the identity matrix and 0 the null vector. Those eigenvectors form a Jordan chain: By recursion starting from v j1 , we have the components of the eigenvector v jk : We express F(n − 1) in the generalized eigenvector basis: with coecients a jk given by the equation VA = F(n − 1), A formed by the coecients a jk , and V with columns formed by the eigenvectors v jk , which defines a confluent Vandermonde matrix [59]. From the inverse of the confluent Vandermonde matrix [60][61][62], we obtain (11). Thus, using representation (A.8) in equation (A.4) and applying the Jordan chain, we have: whose first component is exactly the expression for F(t): was used to obtain the last line, which can finally be written as (10). If all roots are simple, m j = 1, ∀j, we have the simpler formula: where F s (t) is a vector containing generalized random Fibonacci numbers: 2) and Φ is the companion matrix (A.3) and Θ(t) is: By recursion, we write: (C.4) The generalized random Fibonacci number F s (t) is given by the first component of F s (t), which can be expressed as: The first term is just the deterministic generalized Fibonacci number F (t) (equation (10)) and using results from appendix A on generalized eigenvectors, we develop each factor of the sum in the second term: θ i l η i l (n + w l + s)F (w l − (w l−1 + 1) + (n − i l )) (C.6) and obtain the expression for the generalized random Fibonacci numbers as stated in (16).

(D.2)
We write the first sum as s−n   θ i l θ i l η i l η i l F (w l + n − i l )F (w l + n − i l ) .

(D.3)
Terms in which w l → ∞ do not contribute to the sum because F (w l → ∞) → 0 (in the unit root case, F (w l → ∞) converges to a constant). We also impose the conv ergence https://doi.org/10.1088/1742-5468/aaf109 J. Stat. Mech. (2019) 013403 of the expression (D. 3), implying that terms where h → ∞ do not contribute either. Therefore, we can write: which can be expressed as in equation (19).
The expression of C(τ ) for RCAR is the same as for AR but amplified by the factor 1 1−γ , from which we derive the MSD (17) and the autocovariance function (18).