How to fit a regression model with ARIMA errors on the seasonally adjusted component of a time series (in R)?

I want to do these two things (combined) with a time series T:

1. forecast the seasonally adjusted component of T (STL used for the decomposition) and "add back" the seasonality (I assume that the seasonal component is unchanging, so I use naïve method for the seasonal component)
2. fit a regression model with ARIMA errors (exogenous regressors included in the formula)

In other words, I want to obtain forecasts using the seasonally adjusted component of T integrating an external predictor and "adding back" the seasonality.

I can do these two operations separately, but I can't get them to work in combination

Here is some toy examples:

library(forecast)
library(tsibble)
library(tibble)
library(tidyverse)
library(fable)
library(feasts)
library(fabletools)

mutate(Time = yearquarter(Time)) %>%
as_tsibble(index = Time)


Example of fit and forecast with seasonally adjusted component of T:

model_def = decomposition_model(STL,
Consumption  ~ season(window = 'periodic') + trend(window = 13),
SNAIVE(season_year),
dcmp_args = list(robust=TRUE))

fit <- us_change %>% model(model_def)

report(fit)

forecast(fit, h=8) %>% autoplot(us_change)


Example of regression model with ARIMA errors (Income as predictor):

model_def = ARIMA(Consumption ~ Income + PDQ(0,0,0))

fit <- us_change %>% model(model_def)

report(fit)

us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income)) forecast(fit, new_data = us_change_future) %>% autoplot(us_change)  These examples work, but I would like to do something like this: model_def = decomposition_model(STL, Consumption ~ season(window = 'periodic') + trend(window = 13), ARIMA(season_adjust ~ Income + PDQ(0,0,0)), SNAIVE(season_year), dcmp_args = list(robust=TRUE)) fit <- us_change %>% model(model_def) report(fit) us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income))

forecast(fit, new_data = us_change_future) %>% autoplot(us_change)


I get this output in the console:

> fit <- us_change %>% model(model_def)
Warning message:
1 error encountered for model_def

>
> report(fit)
Series: Consumption
Model: NULL model
NULL model>


So I tried doing this in decomposition_model:

model_def = decomposition_model(STL,
Consumption  ~ season(window = 'periodic') + trend(window = 13),
SNAIVE(season_year),
dcmp_args = list(robust=TRUE))


No problem with the fit, but now I get an error in the forecast:

> forecast(fit, new_data = us_change_future) %>% autoplot(us_change)
Error in args_recycle(.l) : all(lengths == 1L | lengths == n) is not TRUE