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;
