How to estimate the parameters of the following log-likelihood function?

I would like to estimate the parameters based on the famous Merton model used probability of default modelling:

Suppose firms' logarithmic returns are following the standard normal distribution and return of a normal firm can be written as:

$R_{it}=\sqrt{\rho} F_t + \sqrt{1-\rho}U_{it}$,

where $F$ is a logarithmic return of an economy at time $t$ independent from the firm $i$ and $U$ is a firm-specific return of firm $i$ at time $t$, both $F$ and $U$ are mutually independent and have standard normal distribution $N(0,1)$.

Coefficient $\rho$ represents correlation between the returns of two randomly chosen firms. Assume that default event of firm $i$ at time $t$ occurs, when firm's return falls below some threshold $S$:

$P(Y_{it}=1)=P(R_{it}<S)$,

where Y_{it} is a binary variable indicating default when equal to $1$ and $0$ otherwise. We model $S$ as a function of linear combination of relevant macroeconomic factors $x$

$S=\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt}$

Then we can write probability of default as

$p_{it}=P(\sqrt{\rho} F_t + \sqrt{1-\rho}U_{it}<\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt}) = \psi(\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt})$

where $\psi$ is a standard normal distribution function. If we denote latent factor $f_t$ as an unobserved realization of variable $F_t$, the probability of default conditional on $f_t$ is

$p_i(f_t)=P \left( U_{it}<{{\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt}-\sqrt{\rho}f_t} \over {\sqrt{1-\rho}}} \right) = \psi \left({{\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt}-\sqrt{\rho}f_t} \over {\sqrt{1-\rho}}} \right)$

To simplify the analysis we assunme all borrowers have the same individual probability $p_i$ and all default events are independent, then according to the “law of large numbers” the default rate of the portfolio can be estimated as individual default probability

$P(p(f_t)=p_i(f_t)|F_t=f_t) = 1$

and unconditional probability can by obtained by integration over random factor

$p_t=P(Y_t=1)=\int_{-\infty}^{\infty}p(f_t)\phi(f_t)df_t$,

where $\phi$ is standard normal density function.

The random factor is assumed to be independent between borrowers. The number of defaults $D_t(f_t)$ is assumed to have a binomial distribution, with conditional default probability $p(f_t)$ and the given number of companies $N_t$.

$D(f_t) \approx Bi(N_t,p(f_t))$

The conditional probability of having exactly $d_t$ defaults at time $t$

$P(D_t=d_t|F_t=f_t)= {n_t \choose d_t} p(f_t)^d_t (1-p(f_t))^{n_t-d_t}$

and then unconditional

$P(D_t=d_t)=\int_{-\infty}^{\infty}{n_t \choose d_t} p(f_t)^d_t (1-p(f_t))^{n_t-d_t}\phi(f_t)df_t$

The parameters of the model for conditional probability of default should be estimable by the following log-likelihood function:

$l(\beta, \rho) = \sum_{t=1}^{T} ln \left\{ \int_{-\infty}^{\infty} {n_t \choose d_T} \psi \left( {{\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt}-\sqrt{\rho}f_t} \over {\sqrt{1-\rho}}} \right)^{d_t} \left[ 1-\psi \left( {{\beta_0 + \sum_{j=1}^{N} \beta_j x_{jt}-\sqrt{\rho}f_t} \over {\sqrt{1-\rho}}} \right) \right]^{n_t-d_t} \phi(f_t) df_t \right\}$

Anyone has an idea what kind of model to implement or how to employ the above formula consistently in R? To my mind comes something as probit random (mixed) effects models, but I don't have experience with them yet. Thanks in advance for your help.