Distribution of $ad-bc$

Question

I'm interested in the stochastic process

$$f_t=(ad-bc)_t$$

where $(a,b,c,d)_t$ is governed by the following transition rules:

$$\begin{align}(a,b,c,d) \rightarrow \begin{cases} (a+1,b,c,d) \;\;\;\; \text{ with probability }p_x &\text{ if }ad>bc\\ (a,b,c+1,d) \;\;\;\; \text{ with probability }1-p_x &\text{ if }ad>bc\\ (a,b+1,c,d) \;\;\;\; \text{ with probability }p_y &\text{ if }ad\leq bc\\ (a,b,c,d+1) \;\;\;\; \text{ with probability }1-p_y &\text{ if }ad\leq bc\\ \end{cases} \end{align}$$

with initial state $(1,1,1,1)$.

I want to compute some quantities relating to $f_t$. In particular its distribution at time $t$, $P_t$, and (I think equivalently), the probability at a given time that $f_t$ changes sign. Maybe also $n$-step transition probability densities would be useful for that goal.

Here are some thoughts: The transitions of $f_t$ derived from the above would be $$ f_{t+1} = \begin{cases} f_t + d_t \;\;\;\; \text{ with probability } p_x &\text{ if } f_t > 0\\ f_t - b_t \;\;\;\; \text{ with probability } 1- p_x &\text{ if } f_t>0\\ f_t - c_t \;\;\;\; \text{ with probability } p_y &\text{ if } f_t \leq 0\\ f_t + a_t \;\;\;\; \text{ with probability } 1-p_y &\text{ if } f_t \leq 0\\ \end{cases} $$

An attempt at a recursive relation for its distribution $P_t$ at time $t$: \begin{align} P_t(f_t) &= P(f_{t-1} > 0) \big[ P_{t-1}(f_t-d_{t-1})p_x + P_{t-1}(f_t+b_{t-1})(1-p_x)\big] \\ &+ P(f_{t-1} \leq 0) \big[ P_{t-1}(f_t+c_{t-1})p_y + P_{t-1}(f_t-a_{t-1})(1-p_y)\big] \end{align}

For starters, $P_0(f) = \delta(0)$ since we assume $(a,b,c,d)$ starts at $(1,1,1,1)$. Then, $P_1(1) = (1-p_y)$ and $P_1(-1) = p_y$. Worst case, it's possible to start from $(1,1,1,1)$ and compute all the possible states in $t$ steps, so I know this must be a tractable problem.

Any input would be valued, thanks very much.

Answer 1

Both $\ \big(a_t,b_t,c_t,d_t\big)\ $ and $\ \big(a_t,b_t,c_t,d_t,f_t\big)\ $ are discrete-time, time-homogeneous Markov chains, so you can obtain the distribution of $\ f_t\ $ (at least for moderate values of $\ t\ $) by using the standard procedure for obtaining the state distribution of such chains. In fact, since $\ f_t\ $ is a deterministic function of $\ \big(a_t,b_t,c_t,d_t\big)\ $, once you have the distribution of the latter, you can compute the distribution of $\ f_t\ $ as $$ \mathbb{P}\big(f_t=\varphi\big)=\sum_{\{(\alpha,\beta,\gamma,\delta)\,|\,\alpha\delta-\beta\gamma=\varphi\}}\mathbb{P}\big(\big(a_t,b_t,c_t,d_t\big)=(\alpha,\beta,\gamma,\delta)\big)\ . $$ The only fly in the ointment is that the number of states accessible at time $\ t\ $ from the initial state $\ (1,1,1,1)\ $ at time $\ t=1\ $ is of order $\ t^3\ $. A state with positive probability at time $\ t\ $ must be a quadruple of positive integers that sum to $\ t+3\ $, and there are $\ {t+2\choose3}\ $ such quadruples. Not all such quadruples will have positive probability, however, so the number of accessible states will be somewhat smaller than this. For $\ p_x=0.7\ $, $\ p_y=0.6\ $ and $\ t=100\ $, for instance, I computed the number of accessible states to be $67,548$, whereas $\ {102\choose3}=176,851\ $. I doubt if the number of accessible states is going to depend much, if at all, on the values of $\ p_x\ $ and $\ p_y\ $ as long as neither of them is $0$ or $1$, so $67,548$ is probably a good indication of the number of accessible states for $\ t=100\ $

With a modern computer I don't think you'll have much of a problem carrying out the computation for $\ t\ $ up to a few hundred, but once $\ t\ $ gets into the thousands, you'll probably start running up against the limits of feasible computability.

Let $\ \pi(t,\alpha,\beta,\gamma,\delta)=\mathbb{P}\big(\big(a_{t+1},b_{t+1},c_{t+1},d_{t+1}\big)=(\alpha,\beta,\gamma,\delta)\big)\ $. Then

• if $\ \alpha+\beta+\gamma+\delta\ne t+3\ $,$\ \alpha\le0\ $, $\ \beta\le0\ $, $\ \gamma\le0\ $, or $\ \delta\le0\ $, $$ \pi(t,\alpha,\beta,\gamma,\delta)=0\ ; $$
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $, $\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \beta\gamma+\delta<\alpha\delta\le\beta\gamma+\alpha\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=p_x\pi(t-1,\alpha-1,\beta,\gamma,\delta)\\ &\hspace{1em}+\big(1-p_x\big)\pi(t-1,\alpha,\beta,\gamma-1,\delta)\\ &\hspace{1em}+\big(1-p_y\big)\pi(t-1,\alpha,\beta,\gamma,\delta-1)\ , \end{align} because the three inequalities $ (\alpha-1)\delta>\beta\gamma\ $,$\ \alpha\delta>\beta(\gamma-1)\ $, and $\ \alpha(\delta-1)\le\beta\gamma\ $ will all be satisfied, while the inequality $\ \alpha\delta\le(\beta-1)\gamma\ $ will not;
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $,$\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \max(\beta\gamma+\delta,\beta\gamma+\alpha)<$$\alpha\delta\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=p_x\pi(t-1,\alpha-1,\beta,\gamma,\delta)\\ &\hspace{1em}+\big(1-p_x\big)\pi(t-1,\alpha,\beta,\gamma-1,\delta) \end{align} because the two inequalities $ (\alpha-1)\delta>\beta\gamma\ $ and $\ \alpha\delta>\beta(\gamma-1)\ $ will both be satisfied, while the inequalities $\ \alpha(\delta-1)\le\beta\gamma\ $ and $\ \alpha\delta\le(\beta-1)\gamma\ $ will not;
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $,$\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \beta\gamma+ \alpha< \alpha\delta\le\beta\gamma+\delta\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=\big(1-p_x\big)\pi(t-1,\alpha,\beta,\gamma-1,\delta) \end{align} because the inequality $\ \alpha\delta>\beta(\gamma-1)\ $ will be satisfied, while the inequalities $ (\alpha-1)\delta>\beta\gamma\ $ , $\ \alpha(\delta-1)\le\beta\gamma\ $ and $\ \alpha\delta\le(\beta-1)\gamma\ $ will not;
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $,$\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \beta\gamma-\beta<\alpha\delta\le\beta\gamma-\gamma\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=\big(1-p_x\big)\pi(t-1,\alpha,\beta,\gamma-1,\delta)\\ &\hspace{1em}+\big(1-p_y\big)\pi(t-1,\alpha,\beta,\gamma,\delta-1)\\ &\hspace{1em}+p_y\pi(t-1,\alpha,\beta-1,\gamma,\delta)\ , \end{align} because the three inequalities $\ \alpha\delta\le(\beta-1)\gamma\ $,$\ \alpha\delta>\beta(\gamma-1)\ $, and $\ \alpha(\delta-1)\le\beta\gamma\ $ will all be satisfied, while the inequality $ (\alpha-1)\delta>\beta\gamma\ $ will not;
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $,$\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \alpha\delta\le\min(\beta\gamma-\gamma,\beta\gamma-\beta)\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=\big(1-p_y\big)\pi(t-1,\alpha,\beta,\gamma,\delta-1)\\ &\hspace{1em}+p_y\pi(t-1,\alpha,\beta-1,\gamma,\delta)\ , \end{align} because the two inequalities $\ \alpha\delta\le(\beta-1)\gamma\ $ and $\ \alpha(\delta-1)\le\beta\gamma\ $ will both be satisfied, while the inequalities $\ \alpha\delta>\beta(\gamma-1)\ $ $ (\alpha-1)\delta>\beta\gamma\ $ will not;
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $,$\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \beta\gamma-\gamma< \alpha\delta\le\beta\gamma-\beta\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=\big(1-p_y\big)\pi(t-1,\alpha,\beta,\gamma,\delta-1) \end{align} because the inequality $\ \alpha(\delta-1)\le\beta\gamma\ $ will be satisfied, while the inequalities $ (\alpha-1)\delta>\beta\gamma\ $ , $\ \alpha\delta>\beta(\gamma-1)\ $ and $\ \alpha\delta\le(\beta-1)\gamma\ $ will not;
• if $\ \alpha+\beta+\gamma+\delta= t+3\ $,$\ \alpha\ge1\ $, $\ \beta\ge1\ $, $\ \gamma\ge1\ $, $\ \delta\ge1\ $, and $\ \max(\beta\gamma-\beta,\beta\gamma-\gamma)<$$\alpha\delta\le$$\min(\beta\gamma+\delta,\beta\gamma+\alpha)\ $, \begin{align} \pi(t,\alpha,\beta,\gamma,\delta)&=\big(1-p_x\big)\pi(t-1,\alpha,\beta,\gamma-1,\delta)\\ &\hspace{1em}+\big(1-p_y\big)\pi(t-1,\alpha,\beta,\gamma,\delta-1) \end{align} because the two inequalites $ \alpha\delta>\beta(\gamma-1)\ $ and $ \alpha(\delta-1)\le\beta\gamma\ $ will both be satisfied, while the inequalites $\ (\alpha-1)\delta>\beta\gamma\ $ and $\ \alpha\delta\le(\beta-1)\gamma\ $ will not.

In practice, computing $\ \pi(s,\cdot,\cdot,\cdot,\cdot)\ $ successively for $\ s=1,2,\dots,t\ $ is actually somewhat simpler than the above tabulation of eight separate cases might appear to suggest.

I wrote the script below for the online Magma calculator. If you copy and paste it into that calculator and press the submit button, it should print out the distribution of $\ f_{20}\ $ for the chosen values of $\ p_x\ $ (viz. $0.7$) and $\ p_y\ $ (viz. $0.6$).

  t:=20;  
  px:=0.7;  
  py:=0.6;  
  I4:=CartesianPower(Integers(),4);  
  newstatedist:=AssociativeArray(I4);  
  oldstatedist:=AssociativeArray(I4);  
  oldstatedist[<1,1,1,1>]:=1.0;  
  for s in [2..t] do  
     for k-> v in oldstatedist do  
        if k[1]*k[2] gt k[3]*k[4] then
          if IsDefined(newstatedist,<k[1]+1,k[2],k[3],k[4]>) then
              newstatedist[<k[1]+1,k[2],k[3],k[4]>]:=newstatedist[<k[1]+1,k[2],k[3],k[4]>]+px*v;
           else
              newstatedist[<k[1]+1,k[2],k[3],k[4]>]:=px*v;
           end if;
           if IsDefined(newstatedist,<k[1],k[2],k[3]+1,k[4]>) then
             newstatedist[<k[1],k[2],k[3]+1,k[4]>]:=newstatedist[<k[1],k[2],k[3]+1,k[4]>]+(1-px)*v;
           else
              newstatedist[<k[1],k[2],k[3]+1,k[4]>]:=(1-px)*v;
           end if;
         else
           if IsDefined(newstatedist,<k[1],k[2]+1,k[3],k[4]>) then
             newstatedist[<k[1],k[2]+1,k[3],k[4]>]:=newstatedist[<k[1],k[2]+1,k[3],k[4]>]+py*v;
           else
             newstatedist[<k[1],k[2]+1,k[3],k[4]>]:=py*v;
           end if;
           if IsDefined(newstatedist,<k[1],k[2],k[3],k[4]+1>) then
              newstatedist[<k[1],k[2],k[3],k[4]+1>]:=newstatedist[<k[1],k[2],k[3],k[4]+1>]+(1-py)*v;
           else
             newstatedist[<k[1],k[2],k[3],k[4]+1>]:=(1-py)*v;
           end if;
         end if;
       end for;  
       oldstatedist:=newstatedist;  
       for k-> v in newstatedist do  
          Remove(~newstatedist,k);  
       end for;  
   end for;  
   f:=AssociativeArray(Integers());  
   for k->v in oldstatedist do  
     if IsDefined(f,k[1]*k[2]-k[3]*k[4]) then  
       f[k[1]*k[2]-k[3]*k[4]]:=f[k[1]*k[2]-k[3]*k[4]]+v;  
     else  
       f[k[1]*k[2]-k[3]*k[4]]:=v;  
     end if;  
    end for;  
    for i->v in f do  
      print i,v;
    end for;

Answer 2

Although not a full solution, here are some thoughts that may help along the way:

An elementary approach would be to formalize the problem and introduce the probability space $$ \Omega = (\mathbb{N}^4)^\mathbb{N}, $$ i.e. the space of sequences $(a_t, b_t, c_t, d_t)_{t\in\mathbb{N}}$, so that the transition rules define a probability measure $\mathbb{P}$ on the power set $\mathcal{P}(\Omega)$.

Then each $f_t$, $t\in\mathbb{N}$, is a random variable from $\Omega \to \mathbb{Z}$ (equivalently, $f$ is a stochastic process from $\Omega \times \mathbb{N} \to \mathbb{Z}$).

Regarding the distribution of $f_t$ we may then write, for every $k\in\mathbb{Z}$, \begin{align} \mathbb{P}(f_t = k) = \mathbb{P}(&\{f_{t-1} = k - d_{t-1} \text{ and } a_t = a_{t-1}+1 \text{ and } f_{t-1} > 0\}\\ &\cup \{f_{t-1} = k + b_{t-1} \text{ and } c_t = c_{t-1}+1 \text{ and } f_{t-1} > 0\}\\ &\cup \{f_{t-1} = k + c_{t-1} \text{ and } b_t = b_{t-1}+1 \text{ and } f_{t-1} \leq 0\}\\ &\cup \{f_{t-1} = k - a_{t-1} \text{ and } d_t = d_{t-1}+1 \text{ and } f_{t-1} \leq 0\}). \end{align} The unions above are clearly disjoint.