Lyapunov exponent and variance in the CLT for products of random matrices related to random Fibonacci sequences

We consider three matrix models of order 2 with one random entry $\epsilon$ and the other three entries being deterministic. In the first model, we let $\epsilon\sim\textrm{Bernoulli}\left(\frac{1}{2}\right)$. For this model we develop a new technique to obtain estimates for the top Lyapunov exponent in terms of a multi-level recursion involving Fibonacci-like sequences. This in turn gives a new characterization for the Lyapunov exponent in terms of these sequences. In the second model, we give similar estimates when $\epsilon\sim\textrm{Bernoulli}\left(p\right)$ and $p\in [0,1]$ is a parameter. Both of these models are related to random Fibonacci sequences. In the last model, we compute the Lyapunov exponent exactly when the random entry is replaced with $\xi\epsilon$ where $\epsilon$ is a standard Cauchy random variable and $\xi$ is a real parameter. We then use Monte Carlo simulations to approximate the variance in the CLT for both parameter models.


Introduction
The main purpose of our paper is to develop new methods to obtain precise estimates of Lyapunov exponents and the variance for the CLT related to the products of random matrices. Let {Y i } i≥1 be a sequence of i.i.d. random matrices distributed according to a probability measure µ. Further, let S n = Y n Y n−1 · · · Y 2 Y 1 .
Assuming that E log + Y 1 < ∞, the top Lyapunov exponent λ associated with µ is given by (1) λ := lim n→∞ 1 n E log S n with λ ∈ R ∪ {−∞}. The top Lyapunov exponent gives the rate of exponential growth of the matrix norm of S n as n → ∞. Since all finite-dimensional norms are equivalent, λ is independent of the choice of norm · . Although λ depends on µ, we usually omit this dependence from our notation. While one can also define a spectrum of Lyapunov exponents, in this paper we will only be concerned with the top Lyapunov exponent λ and we refer to it as simply the Lyapunov exponent. Occasionally, when we are considering λ over a family of distributions parametrized by some variable, we will write λ as a function of that variable. Furstenberg and Kesten (1960) and Le Page (1982) found analogues of the Law of Large Numbers and Central Limit Theorem, respectively, for the norm of these partial products. Despite these results having been established for some time, in most cases it is still impossible to compute the Lyapunov exponent explicitly from the distribution of the matrices. Moreover, computing the variance in the CLT has received scant attention in the literature. We point out that because of the difficulty in computing Lyapunov exponents, most authors need to develop new techniques for specific matrix models rather than work in a general framework.
In this paper, we investigate the behavior of the Lyapunov exponent as the common distribution of the sequence of random matrices varies with a parameter. While there are works in the literature where explicit expressions have been obtained for some matrix models under certain conditions [4,5,6,17,18,19], besides a few special examples, it is not possible to find a general explicit formula for the Lyapunov exponent. There is, however, an extensive literature on approximating the Lyapunov exponent for models where it cannot be calculated explicitly (see [23,22]). For instance, in [22], λ is expressed in terms of associated complex functions and a more general algorithm to numerically approximate λ is given. The method is efficient and converges very fast. The method also applies to a large class of matrix models. There is also a significant interest in computing Lyapunov exponents in physics, with some recent work found in [1,2,7,8,15,16]. The analytic properties of the Lyapunov exponent as a function of the transition probabilities are studied in [20,21,24]. Lyapunov exponents are also useful in mathematical biology in the study of population dynamics.
A random Fibonacci sequence g 0 , g 1 , g 2 . . . is defined by g 0 = g 1 = 1 along with the recursive relation g n+1 = g n ± g n−1 (linear case) or g n+1 = |g n ± g n−1 | (non-linear case) for all n ∈ N, where the sign ± is chosen by tossing a fair or biased coin (positive sign has probability p). In [25], Viswanath studied the exponential growth of |g n | as n → ∞ in the linear case with p = 1 2 by connecting it to a product of random matrices and then employing a new computational method to calculate the Lyapunov exponent to any degree of accuracy. The method involves using Stern-Brocot sequences, Furstenberg's Theorem (see Theorem 2.3) and the invariant measure to compute λ. We also point to the work of [14,13,12] where the authors generalized the results of Viswanath by letting 0 < p ≤ 1 and treating λ as a function of p which bears some similarity to the model we study in Section 3. They also considered the non-linear case.
The model that is most relevant to our results is given in [11], where the authors give an explicit formula for the cumulative distribution function of a random variable X p on (0, ∞) characterized by the distributional identity where p is a Bernoulli (p) random variable independent of X p . Let CDF denote the cumulative distribution function for a random variable. The CDF of X p is given in terms of a continued fraction expansion. We will later see that the distribution of X p is the invariant distribution for the product of random matrices studied in Section 3. We summarize the main results of the paper as follows. Consider the random matrices (2) Lyapunov exponent when ∼ Bernoulli 1 2 (See Theorem 3.2): The Lyapunov exponent λ can be estimated by p n ≤ λ ≤ q n , where p n = log c n (n + 7) 2 n and q n = log c n (n + 4) 2 n , and c n is given by Definitions 3.1 and 3.2. Moreover, lim n→∞ p n = lim n→∞ q n = λ.
The method we develop differs from that of the papers listed above and requires the study of an interesting multi-level recursion satisfied by c n .
(3) Exact Lyapunov exponent involving Cauchy random variable (See Proposition 4.1): When then the Lyapunov exponent λ(ξ) is given by (4) Variance Simulation (See Figures 4, 5 and 6) The paper is organized as follows. In Section 2 we give the preliminaries needed for the paper. In Section 3, we provide exact upper and lower bounds on the Lyapunov exponent associated with the product of random matrices where one entry is Bernoulli (p) with 0 < p < 1. In particular, in Section 3.3 we study the p = 1 2 case and provide a sequence of progressively better bounds. We prove that these bounds converge to the Lyapunov exponent which gives a new characterization for the Lyapunov exponent. Not surprisingly, these bounds are related to Fibonacci sequences as in the work of [11,14,13,12,25].
In Section 4, we give an example of a well-known model where we can calculate the Lyapunov exponent explicitly. In this model, one entry in the random matrix has the Cauchy distribution. In Section 5, we examine the less studied variance associated with a multiplicative Central Limit Theorem for products of random matrices. The multiplicative CLT holds under some reasonable assumptions, see [4]. It states that for converge weakly to a Gaussian random variable with mean 0 and variance σ 2 > 0 as n → ∞. In the special case where the distribution of Y 1 x / x doesn't depend on x ∈ R d \ {0}, Cohen and Newman [6] gave the explicit formulas that hold whenever the expectations are finite. As far as the authors know, this is the only case where an explicit formula for the variance is given. Compared to the calculation of the Lyapunov exponent, there have been relatively few attempts to explicitly compute or numerically approximate the variance. We address this deficiency in the context of the parameter models that we consider by first describing an easy to implement Monte Carlo simulation scheme and then using it to approximate the variance for some of the models we considered earlier in the paper.

Preliminaries
In what follows, we introduce notational conventions and terminology and recall well-known results regarding the Lyapunov exponent. Let P 1 (R) denote the one-dimensional projective space. Recall that we can regard P 1 (R) as the space of all one dimensional subspaces of R 2 . To describe P 1 (R), let us first define the following equivalence relation ∼ on R 2 \ {0}. We say that the vectors x, x ∈ R 2 \ {0} are equivalent, denoted by x ∼ x , if there exists a nonzero real number c such that x = cx . We definex to be the equivalence class of a vector x ∈ R 2 \ {0}. Now we can define P 1 (R) as the set of all such equivalence classesx. We can also define a bijective map φ : Hence with a slight abuse of notation we can Consider the following group action of GL(2, R) on P 1 (R). For A = a b c d ∈ GL(2, R) and x ∈ P 1 (R), we define Let µ and ν be probability measures on GL(2, R) and P 1 (R), respectively. We say that ν is µ-invariant if it satisfies for all bounded measurable functions f : P 1 (R) → R. Furthermore, we say that a set G ⊂ GL(2, R) is strongly irreducible if there is no finite family V 1 , . . . , V k of proper 1-dimensional vector subspaces of R 2 such that For a real valued function f , define f + = max {f, 0}. The following result by Furstenberg and Kesten in [9] gives an important analogue to the Law of Large Numbers. For the rest of this paper, we will suppose that µ is a probability measure on the group GL(2, R) and that the matrices {Y i } i≥1 are distributed according to µ. However, Theorems 2.2, 2.3 and 2.4 all have statements valid for matrices in GL(d, R) as well. In [10], Furstenberg and Kifer give an expression for λ in terms of µ and the µ-invariant probability measures ν on P 1 (R). The following result is given in [10, Theorem 2.2].

Theorem 2.2 (Furstenberg-Kifer)
Let µ be a probability measure on the group GL(2, R) and {Y i } i≥1 be a sequence of i.i.d. random matrices distributed according to µ. If E log + Y 1 + log + Y −1 1 < ∞, then the Lyapunov exponent is given by where the supremum is taken over all probability measures ν on P 1 (R) that are µ-invariant.
If ν is the unique µ-invariant probability measure on P 1 (R), then Theorem 2.2 implies that the Lyapunov exponent can be written as Sufficient conditions for the existence of such a unique ν were given by Furstenberg and can be found in [4, Theorem II.4.1].

Theorem 2.3 (Furstenberg)
Let µ be a probability measure on the group GL(2, R) and {Y i } i≥1 be a sequence of i.i.d. random matrices distributed according to µ. Additionally, let G µ be the smallest closed subgroup containing the support of µ. Suppose the following hold: Then there exists a unique µ-invariant probability measure ν on P 1 (R) and λ > 0. Moreover, ν is atomless. Consequently, where ν is the unique µ-invariant probability measure on P 1 (R). Hence, if X is a random variable distributed according to ν, then Moreover, if A and X are independent, we can also conclude that A · X has the same distribution as X, which we write as A · X ∼ X. This follows from the definition of µ-invariance. Thus, a random variable X with law given by the unique µ-invariant distribution on P 1 (R) must satisfy where a and X are independent. Likewise, the law of any P 1 (R)-valued random variable X which satisfies (5) is µ-invariant hence it must be ν. We make use of this distributional identity for the µ-invariant distribution in later sections.
The following result by Le Page can be found in [4, Theorem V.5.4] and gives a less-studied analogue to the Central Limit Theorem.
. Let µ be a probability measure on the group GL(2, R) and {Y i } i≥1 be a sequence of i.i.d. random matrices distributed according to µ. Moreover, let G µ be the smallest closed subgroup containing the support of µ. Suppose the following hold: Then there exists σ > 0 such that for any converge weakly as n → ∞ to a Gaussian random variable with mean 0 and variance σ 2 .
We remark that the relatively recent paper [3] has relaxed the exponential moment condition (i) to a second moment condition which cannot be improved. In Section 5, we use Monte Carlo simulations to approximate the value of σ 2 for two matrix models that satisfy the hypotheses of Theorem 2.4.

Bernoulli (p) Parameter Model
In this section we consider a random matrix model where the random entry follows a Bernoulli (p) distribution and the parameter of interest is p. Recall that a random variable ∼ Bernoulli (p) if P ( = 1) = p and P ( = 0) = 1 − p. Let µ p be the probability measure on GL(2, R) given by It is straightforward to verify that µ p satisfies hypotheses (i)-(iv) of Theorem 2.3. We verify them here for completeness. For (i), we see that E log + Y 1 < ∞ since p has finite support. For (ii), consider the subgroup G generated by the possible realizations of (6). Since the determinant of each realization has absolute value 1, so to does every matrix in G. Clearly, the closure of G, call itḠ, is a closed subgroup that contains the support of µ p . Hence G µp ⊂Ḡ. Moreover, since the absolute value of the determinant is continuous, every matrix inḠ also has determinant with absolute value 1. It follows that the same holds for G µp as required.
For (iii), we first let F 0 , F 1 , F 2 , F 3 , . . . be the usual Fibonacci sequence 0, 1, 1, 2, 3, 5, . . . Then a simple calculation shows that for each positive integer n, we have Since the powers of the matrix (6) with p = 1 must be in G µp and the norm of the powers grow arbitrarily large with large n, it follows that G µp is unbounded and hence not compact. Lastly, hypothesis (iv) can be checked by way of an equivalent condition given in [4,Proposition II.4.3]. This condition is met as long as for anyx ∈ P 1 (R), the set Sx = M ·x : M ∈ G µp has more than two elements. To see that this holds, suppose at least one of x 1 , x 2 ∈ R is nonzero and consider Drawing the matrix M from (6) with p = 1, we have Since for any x, each of these elements in Sx is distinct, it follows that hypothesis (iv) holds. Since µ p satisfies hypotheses (i)-(iv) of Theorem 2.3, we know there exists a unique µ p -invariant distribution ν p that satisfies (3) and that ν p is atomless. Then by (5), any random variable X p with law ν p must satisfy the distributional identity where p ∼ Bernoulli (p) and is independent of X p . Likewise, the law of any P 1 (R)-valued random variable X p which satisfies (7) is µ p -invariant hence it must be ν p . Using (7) and the fact that ν p is atomless, it is not hard to see that X p ∈ (0, ∞) almost surely. See Goswami [11] for this fact and other facts about X, including an expression for its cumulative distribution function in terms of a continued fraction expansion. In Figures  1A and 1B we show the empirical distribution of 100 000 independent draws from ν 1/2 and remark that the fractal nature of this probability measure is clearly apparent. Let λ(p) be the Lyapunov exponent related to µ p . Using (4) and the fact that X p is non-negative, we can write the Lyapunov exponent associated with µ as 3.1. The general 0 < p < 1 case. In this subsection we study λ(p) for general 0 < p < 1 and obtain two sided bounds depending on the parameter p. First we prove some identities for E [log X p ]. We begin by establishing an identity for E [log X p ] which will be later generalized for the p = 1 2 case and used in proving a limiting result.
Lemma 3.1 If X p is a P 1 (R)-valued random variable satisfying (7), then Proof. Let X p be a random variable satisfying (7). Then X p has law ν p given by Theorem 2.3 applied to random matrices of the form (6). Consequently, we have that 0 < λ(p) < ∞ and it follows from (8) that E [log X p ] is positive and finite. Using (7), we start by writing Adding E [log X p ] to both sides of (9) and dividing by 2 results in Continuing in a similar fashion with (10), we obtain where we use (10) in the last equality. Subtracting 1 − 3p 2 E [log X] from both sides of (11) leads to completing the proof.
Lemma 3.2 If X p is a P 1 (R)-valued random variable satisfying (7), then Proof. Recalling that the distribution of X p has non-negative support, observe that This proves (12) which, along with the fact that the distribution of X p is atomless, allows us to write which proves (13). Combining these two identities now leads to (14).
Next we use these results to establish bounds on the Lyapunov exponent which are dependent on p.
Theorem 3.1 Let µ p be the probability measure on GL(2, R) given by (6). Then the Lyapunov exponent λ(p) associated with µ p can be estimated by Proof. Beginning with the upper estimate, first note that log(2x + 1) ≤ log(3x) for x ≥ 1. Now using Lemma 3.1 and (13), we can write Subtracting 1 3 E [log X p ] from both sides of (15) and recalling (8) leads to the desired result. For the lower estimate, we proceed similarly, noting that log(2x + 1) ≥ log(3x) for 0 < x ≤ 1 and using (14) instead of (13) to write Now the lower bound follows from a simple rearrangement of (16).

Approximating λ(p) by simulation.
Let {Y i } i≥1 be an i.i.d. sequence drawn from µ p , and for some x ∈ R 2 with x = 1, construct Hence it follows from Theorem 2.1 that we can approximate λ by the right-hand side of (17) with n large. Since the log U i terms aren't growing with i, this avoids numerical overflow issues and makes for a robust Monte Carlo scheme.
In Figure 2, we plot simulations for λ(p) in black and the upper and lower bounds from Theorem 3.1 in blue. We discretize [0, 1] into sub-intervals of length 0.01 and use n = 1 000 000 in the Monte Carlo scheme described above. 3.3. The p = 1 2 case. In this section we study λ := λ 1 2 in more detail. To set notation, recall that a random variable ∼ Bernoulli 1 2 if P ( = 1) = P ( = 0) = 1 2 . The probability measure µ on GL(2, R) that we consider is given by We know by the general p case that there exists a unique µ-invariant distribution ν that satisfies (3) and that ν is atomless. Then by (5), any random variable X with law ν must satisfy the distributional identity where ∼ Bernoulli 1 2 and is independent of X. Using (4) and the fact that X is non-negative, we can write the Lyapunov exponent associated with µ as (20) λ = E [log X] .
Unlike in the general case, we will be able to obtain a sequence of upper and lower bounds that converge to λ. Recall that by Lemma  . . .
The string of identities above is obtained by iteratively exploiting the distributional equivalence of X and 1 X + , the independence of X and , and elementary logarithmic identities. We will later see that an interesting pattern emerges. At the first step of the iteration, we are looking at the expected value of the log of one affine function of X that is obtained by taking the inner product of the vector (2, 1) and the vector (X, 1). As we move to the second step of the iteration, we encounter the expectation of the log of the product of two affine functions of X. The first one is obtained by taking the inner product of (3, 2) and (X, 1), while the second is obtained by taking the inner product of (1, 2) and (X, 1). At the third step, we encounter the expected value of the log of the product of four = 2 3−1 affine functions of X; these are obtained by respectively taking the inner product of (X, 1) with the vectors (5, 3), (3, 1), (2,3), and (2,1).
In what follows, we represent the vectors generating the aforesaid affine functions of X via inner products with (X, 1), which we call "coefficient pairs", in an array where the row number corresponding to the n th step of the iteration is n − 1. The first four rows of the array are shown below. We use the symbol → to map the collection of coefficient pairs to the real number representing the product of the sum of entries in each coefficient pair in the row; we make extensive use of these quantities later on. For the k th coefficient pair in row n, let a k n denote the first element and b k n the second. To illustrate this notational convention, consider the example 1 14 E [log (3X + 2) (X + 2)] from (22). This is in row n = 1, so we would refer to the 3 in (3X + 2) as a 1 1 and the 2 as b 1 1 . Similarly, the coefficient of X in (X + 2) would be labeled a 2 1 and the 2 would be labeled b 2 1 . In terms of a k n and b k n , the expression is 1 14 E log a 1 1 X + b 1 1 a 2 1 X + b 2 1 . Now we can define the multi-level recursion that describes the array given in (23). We observe several conspicuous patterns in (23) which are implicit in Definition 3.1. For instance, row n is made up of 2 n pairs and the second half of row n is simply row n − 1 where the elements within the coefficient pairs have been switched. One property that will prove useful is the fact that the first coefficient pair in each row dominates the other pairs occurring in that row in the sense that (24) a 1 n ≥ a k n and b 1 n ≥ b k n for all 1 ≤ k ≤ 2 n . This follows from the recursion in Definition 3.1 and induction on n.
To exhibit a less obvious pattern, we first recall that a "Fibonacci-like sequence" of numbers f 0 , f 1 , f 2 . . . is a sequence determined by the initial values f 0 , f 1 such that f n+1 = f n + f n−1 for all n ∈ N. When f 0 = 0, f 1 = 1, we recover the standard Fibonacci sequence. Fibonacci-like sequences can be expressed by an explicit formula. Let f n (f 0 , f 1 ) represent the nth term in the sequence given initial values f 0 , f 1 . If Now note that given n ∈ N and k ∈ 1, . . . , 2 n−1 , we have a k n+1 = a k n + b k n = a k n + a k n−1 Thus, for each k, the sequences {a k n } and {b k n } will be Fibonacci-like sequences in n for n large enough. We use these observations to help establish bounds on the Lyapunov exponent. In order to find suitable estimates, we first need to establish some preliminary results. These involve proving the string of identities given in (22). We also need to prove some elementary inequalities involving the logarithm of the polynomials given inside the expectations in (22).
First, we extend the identities given in (22) to all n. Proof. We begin with n = 0. By Lemma 3.1 with p = 1 2 we have, Now suppose (26) holds for n. We shall show that (26) holds for n + 1. Note that Moving the last term on the right-hand side of (27) to the left leads to Here we have combined and simplified the products appearing in (27) by using the recursion from Definition 3.1. The result now follows by induction.
We now prove the elementary inequalities needed to estimate (26).
Conversely, when 0 < x ≤ 1, Proof. Note that when x ≥ 1, we have a k n x + b k n ≤ x(a k n + b k n ). Taking products and the log of both sides gives us the desired result. The proof of the 0 < x ≤ 1 case follows similarly.
Using (28) and (29), we can prove that the Lyapunov exponent is bounded by terms dependent only on n. First, we define the following quantities that appear as the rightmost entries of (23).
Definition 3.2 For each n ∈ Z ≥0 , let c n be the product of the sums of coefficient pairs in row n of (23). That is, For example, c 0 , . . . , c 3 are displayed in (23). We remark that the recursion from Definition 3.1 implies (30) c n = c n−1 Now we can state our main result of this section.
Theorem 3.2 Let µ be the probability measure on GL(2, R) given by (18). Then for each n ∈ Z ≥0 , the Lyapunov exponent λ associated with µ can be estimated by where (32) p n = log c n (n + 7) 2 n and q n = log c n (n + 4) 2 n .
Subtracting the last term on the right-hand side of (33) from both sides while recalling (20) leads to For the lower bound, we can repeat this same procedure using (29) and (14) instead of (28) and (13) to arrive at Similarly, this implies log c n (n + 7)2 n ≤ λ.
We now show that these bounds converge to the Lyapunov exponent as n → ∞. We first point out the crude estimate c n ≤ (F n+4 ) 2 n where {F n } := {f n (0, 1)} is the usual Fibonacci sequence. This follows from (30), (24), and the fact that a 1 n = F n+3 for all n ≥ 0. Also note that the well-known asymptotic We end this section with the following two remarks.
Remark 3.1. There doesn't seem to be an obvious recursion among the c n values. In order to compute c n using its definition, we must consider 2 n coefficient pairs. We are able to compute p 25 ≈ 0.204266 and q 25 ≈ 0.225397 but going beyond n = 25 exceeds our computing power. After implementing a simple numerical scheme to compute E [log X] using the CDF of X from Theorem 5.2 of [11] along with (14), we expect that λ ≈ 0.2165.
Remark 3.2. The bounds in Lemma 3.1 from the general p case are analogous to p 0 and q 0 from (32) of the Bernoulli 1 2 model. While we can attempt to improve these bounds by mimicking the proof of Theorem 3.2, unlike in that case, there doesn't appear to be a nice expression for the corresponding bounds p n and q n as n gets larger.

ξ · Cauchy Parameter Model
The parameter model studied in this section is based on the standard Cauchy distribution (that is, Cauchy with location x 0 = 0 and scale γ = 1). Recall that the probability density function of a Cauchy (x 0 , γ) random variable with location x 0 ∈ R and scale γ > 0 is Let µ ξ be the probability measure on GL(2, R) given by The fact that µ ξ satisfies the hypotheses of Theorem 2.3 can be seen through a similar analysis as done in the beginning of Section 3 with some slight differences which we now point out. To verify hypothesis (i), we can use the Frobenius matrix norm to arrive at E[log is the density for Cauchy (0, 1). By elementary computations, this integral is seen to be finite for all ξ. Hypothesis (ii) can be verified in the same manner as for the Bernoulli (p) model. Hypothesis (iii) follows from the unbounded support of . For hypothesis (iv), we can again use the equivalent condition given in [4,Proposition II.4.3]. More specifically, draw M from (35) with = 1 ξ and proceed as in the beginning of Section 3. Hence we know there exists a unique µ ξ -invariant distribution ν ξ such that a random variable X ξ has law ν ξ if and only if it satisfies the distributional identity where ∼ Cauchy (0, 1) and is independent of X ξ . The goal of this section is to find the explicit value of the Lyapunov exponent λ(ξ) related to µ ξ . Following the method from [4, pp. 35], we have an explicit formula for the Lyapunov exponent in terms of the parameter ξ. This will allow us to to study the variance in the Central Limit Theorem related to the products of random matrices of the form (35) as formulated in Theorem 2.4. Since the Lyapunov exponent used in our Monte Carlo simulation scheme will be exact, we can obtain a better approximation for the variance compared to the other models we study.
Proposition 4.1 Let µ ξ be the probability measure on GL(2, R) given by (35). Then the Lyapunov exponent λ(ξ) associated with µ ξ is given by Proof. According to (4), we have λ(ξ) = E log |X ξ | , where X ξ is a random variable satisfying (36). To find the law of such an X ξ , we first guess that it is Cauchy (0, γ) for some γ > 0 and then verify that it satisfies (36) for a particular γ.
Assuming that X ξ ∼ Cauchy (0, γ), the well-known transformation properties of the Cauchy distribution imply that the right-hand side of (36) is also Cauchy distributed, namely Hence ( The proof is complete because of the uniqueness of the distribution v ξ such that (36) is satisfied. Figure 3A shows the graph of λ(ξ) for ξ ∈ [−20, 20]; in Figure 3B, we plot λ(ξ) for ξ ∈ [−1, 1].

Variance Simulation
It is straightforward to verify that the hypotheses of Theorem 2.4 are satisfied for the models we studied in Sections 3 and 4. In fact, much of the reasoning done in the beginning of Sections 3 and 4 to verify the conditions of Theorem 2.3 can be used to verify those of Theorem 2.4. For example, in the Bernoulli (p) model, hypothesis (i) follows from the finite support of µ p . For the Cauchy model, we can again use the Frobenius matrix norm to see that E [exp (t (Y 1 ))] = 2 + ξ 2 x 2 t/2 f (x)dx where f (x) is the density for Cauchy (0, 1). By elementary computations, this integral is seen to be finite when t < 1 and hence hypothesis (i) is also satisfied for this model. Moreover, hypothesis (ii) has already be verified for both models and hypothesis (iii) follows from conditions (ii) and (iii) of Theorem 2.3 which have already been verified.
Thus for 0 < p < 1 and ξ = 0, we know there exists σ(p), σ(ξ) > 0 such that for any x ∈ R 2 \ {0}, 1 √ n log S n x − nλ(p) and 1 √ n log S n x − nλ(ξ) converge weakly as n → ∞ to Gaussian random variables with mean 0 and variance σ 2 (p) and σ 2 (ξ). Here the S n are products of matrices distributed according to the probability measures µ p and µ ξ given in Sections 3 and 4, respectively. Motivated by these considerations and following the idea of Section 3.2, we can approximate σ 2 (p) and σ 2 (ξ) by computing the variance of log U i − nλ(p) and L ξ : with n large. Here, as in Section 3.2, the sequence {U i } i≥0 is constructed recursively by U 0 = x and U i = Y i Ui−1 Ui−1 for some x ∈ R 2 with x = 1 and {Y i } i≥1 an i.i.d. sequence drawn from µ p or µ ξ as appropriate. While we have an exact expression for λ(ξ), we must settle for the approximation of λ(p) obtained by simulation in Section 3.

2.
In what follows, we summarize the simulation procedure for σ 2 (p). The procedure for σ 2 (ξ) is practically identical.
(1) Choose an interval [a, b] as the range of p. Divide this interval into sub-intervals of length k where k divides b − a. Let p be of the form a + jk for j = 0, 1, . . . , b−a k . (2) Choose a unit vector x ∈ R 2 .
(3) Simulate L p for each p from Step 1 and store the result as a data vector of length b−a k + 1. (4) Repeat Step 3 an m number of times to obtain an m × b−a k + 1 matrix, where the j th column contains all of the L p simulations corresponding to p = a + (j − 1)k. (5) Estimate Var L a+(j−1)k by the sample variance of the j th column of the matrix.
Note that in all of our simulations, we set x = in Step 2.
We first approximate the variance for the Bernoulli (p) model considered in Section 3. Trivially, we have that σ 2 (0) = σ 2 (1) = 0. For 0 < p < 1, we simulate Var (L p ) with k = 0.01, n = 1000, and m = 1 000 000. We plot the resulting points in Figure 4 and remark that the graph exhibits distinct asymmetry with the maximum variance occurring around p = 0.56. For the Cauchy parameter model from Section 4, it is clear that σ 2 (0) = 0. For ξ = 0, we simulate Var (L ξ ) over both a large and small range of ξ. Figure 5 illustrates the results for ξ ∈ [−20, 20] with k = 0.25. This is the same interval used to produce Figure 3A. In Figure 6, we plot Var (L ξ ) for ξ ∈ [−1, 1] with k = 0.01 to give a much finer resolution of the graph around the origin.