Passing exogeneous variable to ARIMAX

Hello team,

I have been trying to fit an ARIMAX model with order (0,1,1) in Stan with an exogenous variable. I have the Stan code for ARIMA, and I want to reuse this code to fit ARIMAX passing an exogenous variable to it.

data {
  int<lower=0> T;       //no. of observations
  real y[T];           // observed outputs
parameters {
  real phi;           //AR coeff              
  real theta;          //MA coeff
  real mu;        //intercept coeff           
  real<lower=0> sigma;       // noise scale
model {
  vector[T] nu;              // prediction for time t
  vector[T] err;             // error for time t
  nu[1] = mu + phi * mu;     // assume err[0] == 0
  err[1] = y[1] - nu[1];
  for (t in 2:T) {
    nu[t] = mu +  phi * (y[t-1]) + theta * err[t-1];
    err[t] = y[t] - nu[t];
  mu ~ cauchy(-0.07, 0.1); 
  phi ~ normal(0.85, 0.1);
  theta ~ normal(0, 1);
  sigma ~ cauchy(0, 1);
  err ~ normal(0, sigma);   

Is there any way to pass the exogenous variable to this stan code and if yes, how can I reuse this code to fit the ARIMAX model? Is there some other code that I can use to fit the ARIMAX model with order (0,1,1) in Stan?
All suggestions are welcome.

Exogenous values are supplied either through data or transformed data blocks.
@asael_am’s GitHub - asael697/bayesforecast: Automatic forecasting and Bayesian modeling for time series with Stan seems to have arimax feature, and tracing its implementation might be helpful e.g. this line receives xreg. Another resource is Prophet library, this code or this issue on regressors.