Part 1

#Create time series 

soybeans=tq_get(c("PSOYBUSDM"),get="economic.data",from="2000-01-03")

soybean_ts=ts(soybeans$price,frequency=12,start=c(2000,1))

#Decompose ts 

soybean_decomp=decompose(soybean_ts)

#Plot decomposition

plot(soybean_decomp$trend)

plot(soybean_decomp$seasonal)

plot(soybean_decomp$random)

soybean_decomp_out=soybean_decomp[1:4] %>%

  as_tibble() %>%

  rename(price=x) %>% 

  mutate(measure_date=soybeans$date)

#Extract components 

soybean_trend=soybean_decomp$trend

soybean_seasonal=soybean_decomp$seasonal

soybean_residuals=soybean_decomp$random

#Create seasonal adjusted 

adjstd_soybean=soybean_trend + soybean_residuals

plot(adjstd_soybean)

adjstd_soybean_df=as.data.frame(adjstd_soybean)%>%

  rename(sa_price=x)%>%

  mutate(measure_date=soybeans$date)

#Forecast 3 years 

soybean_trend_forecast=forecast(soybean_trend, h=36)

plot(soybean_trend_forecast)

soybean_seasonal_forecast=forecast(soybean_seasonal,h=36)

plot(soybean_seasonal_forecast)

#Create time series object 

soybean_residuals_ts=ts(soybean_residuals,frequency = 12,start=c(2000,1))

# Fit an MA(1) model to the residuals

soybean_residuals_model <- arima(soybean_residuals_ts, order=c(0,0,1))

#Forecast residual for three years 

soybean_residuals_forecast=forecast(soybean_residuals_model,h=36)$mean

plot(soybean_residuals_forecast)

#Combine three forecasts 

soybean_forecast_combined=soybean_trend_forecast$mean + soybean_seasonal_forecast$mean + soybean_residuals_forecast

plot(soybean_forecast_combined)

#Create data frame 

soybean_forecast_combined_df=tibble(price = soybean_forecast_combined)%>%

  mutate(measure_date=seq(as_date("2023-01-01"),by="months",length.out=nrow(.)))

soybean_forecast_combined_df=soybean_forecast_combined_df %>%

  filter(measure_date > max(soybean_decomp_out$measure_date)) %>%

  mutate(forecast=T)

final_out=bind_rows(soybean_decomp_out,soybean_forecast_combined_df)

write.csv(final_out,"soybean_forecast.csv")

Part 2