Consider a data series that has a seasonal component, a trend, and an ARIMA part. I want to forecast this series.

data_ts <- ts(data, frequency = 24)
data_deseason <- stl(data_ts, t.window=50, s.window='periodic', robust=TRUE)
f <- forecast(data_deseason, method='arima', h = N)

By implementing the above code, I am not able to choose the parameters of the ARIMA part, which I would like to.

I tried splitting the data into a season a trend and a remainder part. But then how do I forecast it?

Should I make an ARIMA model for both the trend and the remainder?

trend_arima <- Arima(data_deseason$time.series[,'trend'], order = c(1,1,1))
remainder_arima <- Arima(data_deseason$time.series[,'remainder'], order = c(1,1,1))

and then use forecast() and add the two above components and the season.

Or is there some way to extract the trend model that STL has found?