<p>Hi stats_noob: I read through your question more carefully so I'm putting some questions in another answer. It's definitely NOT answer to your problem but contains some comments and confusions.</p>
<p>1a) I like the idea of using the predictions ( you call them expected values and they can be thought of this way. I like to think of them as predictions ) to make your decisions. I don't think there are constraints in your problem but I'm not totally sure. If there are, then this type of solution may not be feasible and a dynamic optimization might be a better way to go. This leads to 1b).</p>
<p>1b) Dynamic optimization is used when one has constraints in their problem along with some kind of model. I hope you don't have constraints because it's kind of a macro-econ field in itself and somewhat involved. I can provide some intro material on it if there are constraints. This leads into 2).</p>
<p>2 ) At first I thought that the amount of money that you could spend on fertilizer was constrained and that you would need to spend more when fertilizer cost was high and less when fertilizer cost was low ?
But, after reading the complete question, I don't think that's the case.</p>
<p>But then I asked: &quot;Why is the fertilizer cost changing ?&quot;. I'm not clear on that which kind of makes it difficult to provide any insights. Also, how does fertilizer cost act as a predictor for crop yield ? If the fertilizer cost changes, should one expect that to increase or decrease the crop yield ? what are the units of crop yield and fertilizer cost ? Maybe fertilizer cost is not the <span class="math-container">$F_t$</span> in the model. You used the term &quot;fertilizer&quot; so maybe <span class="math-container">$F_t$</span> denotes the amount of ferilizer used at time <span class="math-container">$t$</span>.</p>
<p>3 ) Also, I'm sure that auto.arima is quite a useful R function for estimation but, for your particular model, it's easier to just make a data.frame out of the 4 predictors and then use that in a call to lm. You're clearly very good with R so I think you can do that without my help. If you want some, let me know. Once you do that, then testing various models is just a straightforward call to lm and then you don't have to worry about all of these arima coefficients lying around. Personally, I still don't know what's under the hood of auto.arima but I know with certainty that you don't need it.</p>
<p>4 ) I realize that your R code is just a template-simulation ( not the real data etc. It's beautiful R code by the way. I've never used tidyverse or ggplot but I can still tell that you know your stuff ) I'm not clear on why you're predicting 13 steps ahead ? Shouldn't you be predicting 1 step ahead, constructing a forecast  and then re-estimating the model for the next period so that you can make another forecast ? I might be confused about this also ?</p>
<p>So, if you can address my issues-questions, maybe I can give some useful feedback instead of confusion !!!!!!!</p>
<p>P.S. : 1 R tip: Don't call your data, &quot;data&quot;. That can be dangerous because there was  a function in Splus called data().
I think there's one in R also so doing that can cause confusion or even more serious problems than that.</p>
<p>#===============================================================
ADDENDUM: ADDING RCODE THAT SHOWS HOW TO USE lm TO ESTIMATE THE MODEL DESCRIBED IN THE ORIGINAL QUESTION. DONT KNOW WHY I HAVE TO SKIP LINES. IF I DO NOT, THINGS GET UGLY.</p>
<pre><code>library(zoo)

set.seed(123)

n &lt;- 10

dates &lt;- seq(as.Date(&quot;2020-01-01&quot;), by = &quot;month&quot;, length.out = n)

fertcost &lt;- runif(n, 1000, 1500)

cropyield &lt;- runif(n, 1000, 1500)

df &lt;- data.frame(dates, fertcost, cropyield)

mydata &lt;- read.zoo(df)

mylaggeddata &lt;- merge(Ytee = mydata[,2], Yteem1 = lag(mydata<span class="math-container">$cropyield, -1,), Yteem2= lag(mydata$</span>cropyield, -2), Fteem1= lag(mydata<span class="math-container">$fertcost, -1,),
                      Fteem2 = lag(mydata$</span>fertcost, -2))

mylaggeddata &lt;- as.data.frame(mylaggeddata)

lm(data = mylaggeddata, as.formula(Ytee ~ Yteem1 + Yteem2 + Fteem1 + Fteem2))
</code></pre>
<p>ML: 08-20-2024
#====================================================================</p>
<p>Hi stats_noob: I hope that the R code helped. Notice that, once the observations are aligned correctly, your model is no longer a time series regression model in the sense that time has been eliminated. Therefore, you just have a regression model which can be estimated using lm. Of course, if the error terms was not white noise and there was a more complicated error term structure, then you couldn't use lm.</p>
<p>I read through your answers to my question and I like the idea of forecasting each month out for 12 months and adding them up. The only problem with that approach is that you need the values of <span class="math-container">$F_{t-1}$</span> and <span class="math-container">$F_{t-2}$</span> which you won't have as you go forward. So, that approach is probably not feasible ?</p>
<p>The other confusion I had was that I'm still not sure if <span class="math-container">$F_t$</span> is amount of fertilizer purchased  or fertilizer cost ? In either case, you need to justify their ability as predictors of yield because I'm not sure of their predictive ability ? Fertilizer cost would make sense if you were spending more for a &quot;more powerful&quot; fertilizer but, other than that, I don't see their predictive ability ?</p>
<p>I hope the R code and my questions help ?</p>

CC BY-SA 4.0
2024-08-19T02:46:45.833Z
@ Mark Leeds: Thank you so much for all your efforts and kind words!
stats_noob
CC BY-SA 4.0
2024-08-19T02:49:13.443Z
About 4) . I am predicting 13 steps ahead because I think this is a realistic way to test models. Farmers usually need to forecast 12 months ahead and allow 12 months worth of uncertainty into their decision making. Of course predict 1 month, re-estimate, predict 1 month, re-estimate would give you much better results... but this would not be realistic in my opinion. In my opinion, a farmer might need to buy 12 months worth of fertilizer at once ... regardless of whether half way through the year the farmer's decision turned out to be sub-optimal
stats_noob
CC BY-SA 4.0
2024-08-19T02:51:51.480Z
it reminds me of a situation in college. one time, a chemistry professor said there is a midterm and a final. exam if you skip the midterm, the final exam will be worth the entire course mark. I innocently told the professor : "on the day of the midterm, can I sit inside, quickly flip through the midterm and see how confident I am .. and then decide if I want to write it or leave"? Incredulous, the prof replied "this is a chemistry midterm, not a buffet! where do you students come up with these ideas? if you spent half the time studying as you do coming up with these ideas you would get an A+"
stats_noob
CC BY-SA 4.0
2024-08-19T02:52:54.383Z
The same way, I imagine that the farmer cant "take a peak at the farm conditions each month and re-update his fertilizer purchases ... return excess amounts and place more fertilizer orders if necessary". The farmer is forced to make a decision for 12 months at once
stats_noob
CC BY-SA 4.0
2024-08-19T02:53:57.507Z
About 2) I wanted to make it that in the past, the farmer had spent varying amounts on fertilizer each month ... but when planning for the future, the farmer decides to split the year's cost over 12 months
stats_noob
CC BY-SA 4.0
2024-08-19T02:54:29.883Z
I would be very curious to hear your lm() approach in R. If you have time, can you please write it?
stats_noob
CC BY-SA 4.0
2024-08-19T02:54:37.850Z
BTW - have you heard of linear state space models?
stats_noob
CC BY-SA 4.0
2024-08-19T09:29:34.680Z
let me print out your comments and get back to you in a day or two. I wrote something in R ( I read your things quickly ) where I do an lm on artificial data. The sad thing is that I don't know how to put it in that window that you put yours in. Can you explain how to do that ? The code is short and simple. You'll see what's going on right away. Thanks.
mark leeds
CC BY-SA 4.0
2024-08-19T09:34:40.487Z
Yes, I have some experience with state space models. But I always try to figure out which variables might be helpful and possible relationships between them BEFORE I decide on a methodology. That's kind of secondary, atleast from  my viewpoint. Others may feel differently. It's all good. But yes, a kalman filter is useful when you don't necessarily believe that your predictors are constant. Basically, the predictors vary over time, possibly. Other than that reason, I see no use for it with your current model.
mark leeds
CC BY-SA 4.0
2024-08-19T09:40:09.907Z
If the window creation thing is complicated, I can just put the R code in another answer ? But that would be my third answer to the same question !!!!! I imagine that's discouraged by whomever are the kind people who run this list ? Shout out to them. It's a wonderful list with many generous and kind people.
mark leeds
CC BY-SA 4.0
2024-08-19T10:04:19.003Z
I managed to make the R code to use lm an addendum in the second answer. Only negative is that I had to take out my comments to get it to look decent. If any of the R code is not clear, just let me know. I have comments in the actual code.
mark leeds
CC BY-SA 4.0
2024-08-20T15:55:19.183Z
@stats_noob: I read through your answers and put another addendum below the addendum containing the R-code. I hope it helps.
mark leeds
CC BY-SA 4.0
2024-08-20T17:15:10.210Z
I have edited the answer. Use curly-brackets symbol '{}' for R-code.
Jesper Hybel
CC BY-SA 4.0
2024-08-21T15:49:39.133Z
Thanks Jesper: I'll go over it to see what you did. It's much appreciated.
mark leeds
CC BY-SA 4.0
2024-08-21T15:57:49.290Z
Jesper: I went into edit mode to see what you did but I can't find the {}. I must have woken up too early this morning. It definitely looks much better so thanks but I can't figure out what you did to do that. uuuugh.
mark leeds
CC BY-SA 4.0
2024-08-19T01:03:50.003Z
2024-08-20T17:14:22.600Z
2024-08-20T17:14:22.600Z
mark leeds
<p>Let me just give a short example as to why two lags of <span class="math-container">$y_t$</span> might be &quot;redundant&quot; if you will.</p>
<p>Suppose you just had an AR(1), so</p>
<p><span class="math-container">$$y_t = \beta y_{t-1} + \epsilon_t$$</span></p>
<p>Then, what is this model REALLY saying. Converting it using the lag operator, one obtains:</p>
<p>(1)
<span class="math-container">$$y_t = \sum_{i=0}^{\infty} \beta^{i} \epsilon_{t-i} $$</span></p>
<p>So, assume a shock of <span class="math-container">$\epsilon_t$</span> at time <span class="math-container">$t$</span> ( this is referred to as the impulse in the econometrics literature ) and assume no shocks after that. Then this model implies that the impulse response of a shock is an exponentially smoothed version of that same shock.  So, let us call that shock <span class="math-container">$\epsilon_t$</span> at time  <span class="math-container">$t$</span> as <span class="math-container">$\epsilon^{*}$</span> and assume that there are no shocks after that.</p>
<p>Then, based on (1), first, at time  <span class="math-container">$t$</span> the effect of the <span class="math-container">$\epsilon^{*}$</span> shock is just <span class="math-container">$\epsilon^{*}$</span>. Then, at time <span class="math-container">$t+1$</span>, the shock's effect get decreased to <span class="math-container">$\beta \epsilon^{*}$</span>. Then at time <span class="math-container">$t+2$</span>, this same shock has effect of <span class="math-container">$\beta^2 \epsilon^{*}$</span> and so on and so forth until the effect of the shock totally dissipates and eventually becomes zero.</p>
<p>Now, suppose the model had been an AR(2), rather than an AR(1). So, in this case, one has:</p>
<p><span class="math-container">$$y_t = \beta y_{t-1} + \theta y_{t-2} + \epsilon_t $$</span></p>
<p>So, going through a similar derivation one obtains <span class="math-container">$$y_t = \sum_{i=0}^{\infty}\theta^{i} \epsilon_{t-i-1}$$</span></p>
<p>So not only do we have this one set of &quot;effects&quot; due to the
shock  <span class="math-container">$\epsilon^{*}$</span> at time <span class="math-container">$t$</span>. We then have this second set of exponentially smoothed effects due to the same <span class="math-container">$\epsilon^{*}$</span>. This set starts at time <span class="math-container">$t+1$</span> and result in sequential effects that are exponentially smoothed over time. The difference between this second set of shocks and the first set is that they start at time <span class="math-container">$t+1$</span> and have a different smoothing coefficient with the coefficient being <span class="math-container">$\theta$</span> rather than <span class="math-container">$\beta$</span>.</p>
<p>My point about lagging once rather than twice is that I find this interpretation for the second set of effects difficult to justify when the shock they originate from is the same shock just lagged by one more time unit. Why would there be this &quot;similar&quot; set of effects except lagged by one time unit ? The first set already takes care of the time lagging effect decreasing over time so there probably shouldn't be the need for a second set ? It seems like too much overlapping to me unless there was some known reason for it.</p>
<p>I provided the argument for the AR(1) but a similar argument could be used for the ARDL that you're using. The <span class="math-container">$F_t$</span> are assumed exogenous so one cannot argue using the <span class="math-container">$\epsilon_t$</span>s like we did in the AR examples above but still a similar argument holds: The effects of the <span class="math-container">$F_t$</span> coefficients also dissipate over time due to the auto-regressive-ness of the model so their set of effects are also exponentially smoothed over time. This is why I would argue against a second lag for the <span class="math-container">$F_t$</span> also.</p>
<p>When I was reading about ARDLs ( I think they are synonymous with what is referred to by statisticians as ARIMAX but I'm not certain ),  I found that there is a lot of good political science literature on distributed lag models. This was somewhat surprising and it helped me a lot. The authors that I highly recommend are Luke Keele, Nathaniel Beck and Suzanne Deboef. If you can't find anything connected to them and distributed lags, let me know and I can send some  links. Some of their papers were real eye-openers for me and even better than most of explanations in the econometrics texts.  All the best.</p>
<p>ADDENDUM TO ORIGINAL POST
#================================================================</p>
<p>Notice that the exponentially smoothed &quot;set&quot; of effects arise
ONLY WHEN THERE IS A LAGGED AUTO-REGRESSIVE term, in this case,
<span class="math-container">$y_{t-1}$</span> and <span class="math-container">$y_{t-2}$</span>. If there were no lagged auto-regressive
terms then, any shock or exogenous effect could be viewed as an impulse but it would be a one time impulse at time <span class="math-container">$t$</span> and then zero after that. So, this &quot;lagged effect&quot; concept from the impulse response would not exist without the existence of atleast one of the lagged auto-regressive terms, <span class="math-container">$y_{t-1}$</span> or <span class="math-container">$y_{t-2}$</span>.</p>
<p>SECOND ADDENDUM TO ORIGINAL POST
#==========================================</p>
<p>Note that when explaining the two &quot;sets of effects&quot; caused by having two lagged versions of the same variable, I simplified things are far as the two sets being seperate and additive. The effect from the first lag is being fed into the second lag so the effects are definitely not simply seperate and additive. One would have use the polynomial equivalence approach used in transfer function modelling to figure out the exact sequence of effects due to the impulse response. Nevertheless, even though the way of thinking about the two sets being &quot;additive&quot; is not correct, the argument still holds: &quot;Why would an impulse response of <span class="math-container">$\epsilon^{*}$</span> cause such a complex relationship for the effects of that response over time. It's possible but not likely.</p>

CC BY-SA 4.0
2024-08-16T21:59:30.637Z
@ Mark Leeds: thank you so much for your detailed answer! Just a question - does the overall idea sound right? we can predict the expected values - and then based on the predicted expected values conduct the cost benefit analysis?
stats_noob
CC BY-SA 4.0
2024-08-17T03:14:01.260Z
Hi: Youre welcome. But that last part is really important as far as understanding that you can just take two seperate effects and look at them individually ( like I did ) because what you really have are two filters in series so  SYSTEM RESPONSE would need to be figured out.  So I simplified things terribly when I made a feeble attempt to explain what's being assumed about the impulse response when an AR(2) is used.
mark leeds
CC BY-SA 4.0
2024-08-17T03:26:55.730Z
Hi stats_noob: Let me read it again over the next couple of days because I really didn't read it the way it should be read. I just focused on the estimation part. So, let me get back to you. It sounds like an important question so I will just say whatever I would do if I was doing what you're doing. Others may have different ways ?  I've realized in the last 10 or so years that the amount of overlap in different fields is enormous but even the overlap varies because it depends on what is emphasized. That's the most interesting thing: There's variation in the overlap itself !!!!!!
mark leeds
CC BY-SA 4.0
2024-08-17T03:32:45.607Z
first sentence of next to last comment is "you can't just take". To figure out the overall response given the AR(2) you have to compute the transfer function. This sort of thing falls more under the topic of "system identification". It is also covered in "signals and systems" type classes. For whatever reason, it's  more emphasized in digital signal processing. But we're getting off track so I'll leave it at that.
mark leeds
CC BY-SA 4.0
2024-08-17T06:39:13.923Z
stats_noob: here's a link to Hyndman's text where he explains how the auto.arima model works. I'm just sending it now so that I don't forget. It would be useful to find an example where auto.arima is used to fit a model like yours rather than a reg model with arima errors. He likes the latter approach which has equivalences to the  ARDL models which are more in line with your notation-setup.  I haven't done this sort of thing ( fit ARDL models ) in a while but the systemfit package might be a better option for estimation of your model.  https://otexts.com/fpp2/dynamic.html
mark leeds
CC BY-SA 4.0
2024-08-16T17:56:07.113Z
2024-08-16T20:30:48.190Z
2024-08-16T20:30:48.190Z
mark leeds
<p>I am working on this economics problem and feel a bit lost. I was hoping someone could help me out.</p>
<ul>
<li>Suppose there is a farm in which there are two variables: Crop Yield and Fertilizer Cost</li>
<li>I am interesting in making a time series model to model the relationship between Crop Yield and Fertilizer Cost (I.e how much Fertilizer I choose to buy, assuming I fix the amount of money I am willing to spend on Fertilizer as constant)</li>
<li>Specifically, I believe that Crop Yield is influenced by the (lag-1) and (lag-2) values of Crop Yield and Fertilizer Cost.</li>
<li>After doing some research, it seems to me that an ARIMAX model is a good choice for this problem.</li>
<li>In the end, I want to answer the following questions:
<ul>
<li><strong>What is the (estimated) impact of each extra dollar of fertilizer cost on crop yield?</strong></li>
<li><strong>Suppose I have 3 options over the next 12 months : Invest 1000 dollars each month on fertilizer, Invest 2000 dollars each month on fertilizer, Invest 3000$ each month on fertilizer. Which of these options will result in the largest (estimated) cumulative crop yield over 12 months?</strong></li>
</ul>
</li>
</ul>
<p>I tried to express this in mathematical terms:</p>
<p>Let <span class="math-container">$ Y_t $</span> represent the crop yield at time <span class="math-container">$ t $</span>, and let <span class="math-container">$ F_t $</span> represent the fertilizer cost at time <span class="math-container">$ t $</span>. We assume that the price of fertilizer remains constant, so variations in <span class="math-container">$ F_t $</span> reflect the quantity of fertilizer purchased.</p>
<p>The crop yield <span class="math-container">$ Y_t $</span> is influenced by:</p>
<ul>
<li>The crop yield at lag-1 (<span class="math-container">$ Y_{t-1} $</span>) and lag-2 (<span class="math-container">$ Y_{t-2} $</span>),</li>
<li>The fertilizer cost at lag-1 (<span class="math-container">$ F_{t-1} $</span>) and lag-2 (<span class="math-container">$ F_{t-2} $</span>).</li>
</ul>
<p>The relationship between crop yield and fertilizer cost can be modeled using an ARIMAX model as follows:</p>
<p><span class="math-container">$$
Y_t = \alpha_0 + \alpha_1 Y_{t-1} + \alpha_2 Y_{t-2} + \beta_1 F_{t-1} + \beta_2 F_{t-2} + \epsilon_t
$$</span></p>
<p>where: <span class="math-container">$ \alpha_0 $</span> is the intercept term. <span class="math-container">$ \alpha_1, \alpha_2 $</span> are the coefficients for the lagged crop yields, <span class="math-container">$ \beta_1, \beta_2 $</span> are the coefficients for the lagged fertilizer costs, <span class="math-container">$ \epsilon_t $</span> is the error term.</p>
<p><strong>Thus, it seems to me that <span class="math-container">$\beta$</span> would represent the estimated impact of fertilizer on crop yield ... and I could add the estimated 12 month forecasts for each option to find out the net benefit? Is my understanding of Economics correct?</strong></p>
<p>Can someone please tell me if I am on the right track?</p>
<p>Thanks!</p>
<ul>
<li>Note: I tried to make an R simulation to represent my problem.</li>
</ul>
<p>First I simulate data:</p>
<pre><code>library(forecast)
library(ggplot2)
library(dplyr)
library(tsibble)

set.seed(123)
n &lt;- 100
dates &lt;- seq(as.Date(&quot;2020-01-01&quot;), by = &quot;month&quot;, length.out = n)
fertilizer &lt;- runif(n, 500, 1500)
seasonal_effect &lt;- sin(2 * pi * (1:n) / 12) * 200
trend &lt;- seq(0, 500, length.out = n)
crop_yield &lt;- 0.5 * fertilizer + seasonal_effect + trend + arima.sim(list(ar = c(0.7, 0.2)), n) + rnorm(n, 0, 50)

data &lt;- tibble(
    date = dates,
    crop_yield = crop_yield,
    fertilizer = fertilizer
) %&gt;%
    as_tsibble(index = date)
</code></pre>
<p>Then, I fit the model:</p>
<pre><code>model &lt;- auto.arima(data<span class="math-container">$crop_yield, xreg = cbind(data$</span>fertilizer))
</code></pre>
<p>From there, I prepared the data for the forecasting:</p>
<pre><code>future_dates &lt;- seq(from = max(dates), by = &quot;month&quot;, length.out = 13)

generate_forecast &lt;- function(cost) {
    future_fertilizer &lt;- c(tail(data$fertilizer, 1), rep(cost, 12))
    forecast_result &lt;- forecast(model, xreg = future_fertilizer, h = 13)
    
    list(
        date = future_dates,
        crop_yield = as.numeric(forecast_result<span class="math-container">$mean),
        fertilizer = future_fertilizer,
        lower = as.numeric(forecast_result$</span>lower[, 2]),
        upper = as.numeric(forecast_result<span class="math-container">$upper[, 2]),
        scenario = paste0("Fertilizer $</span>&quot;, cost)
    )
}
</code></pre>
<p>Then, I fit the forecasts and plotted them:</p>
<pre><code>forecast_1000 &lt;- generate_forecast(1000)
forecast_2000 &lt;- generate_forecast(2000)
forecast_3000 &lt;- generate_forecast(3000)

ggplot() +
    geom_line(data = data, aes(x = date, y = crop_yield, color = &quot;Historical&quot;)) +
    geom_point(data = data, aes(x = date, y = crop_yield, color = &quot;Historical&quot;)) +
    geom_line(data = tibble(date = forecast_1000<span class="math-container">$date, crop_yield = forecast_1000$</span>crop_yield), aes(x = date, y = crop_yield, color = &quot;Fertilizer <span class="math-container">$1000")) +
    geom_line(data = tibble(date = forecast_2000$</span>date, crop_yield = forecast_2000<span class="math-container">$crop_yield), aes(x = date, y = crop_yield, color = "Fertilizer $</span>2000&quot;)) +
    geom_line(data = tibble(date = forecast_3000<span class="math-container">$date, crop_yield = forecast_3000$</span>crop_yield), aes(x = date, y = crop_yield, color = &quot;Fertilizer <span class="math-container">$3000")) +
    geom_ribbon(data = tibble(date = forecast_1000$</span>date, ymin = forecast_1000<span class="math-container">$lower, ymax = forecast_1000$</span>upper), aes(x = date, ymin = ymin, ymax = ymax, fill = &quot;Fertilizer <span class="math-container">$1000"), alpha = 0.2) +
    geom_ribbon(data = tibble(date = forecast_2000$</span>date, ymin = forecast_2000<span class="math-container">$lower, ymax = forecast_2000$</span>upper), aes(x = date, ymin = ymin, ymax = ymax, fill = &quot;Fertilizer <span class="math-container">$2000"), alpha = 0.2) +
    geom_ribbon(data = tibble(date = forecast_3000$</span>date, ymin = forecast_3000<span class="math-container">$lower, ymax = forecast_3000$</span>upper), aes(x = date, ymin = ymin, ymax = ymax, fill = &quot;Fertilizer <span class="math-container">$3000"), alpha = 0.2) +
    scale_color_manual(values = c("Historical" = "black", 
                                  "Fertilizer $</span>1000&quot; = &quot;blue&quot;, 
                                  &quot;Fertilizer <span class="math-container">$2000" = "red", 
                                  "Fertilizer $</span>3000&quot; = &quot;green&quot;)) +
    scale_fill_manual(values = c(&quot;Fertilizer <span class="math-container">$1000" = "blue", 
                                 "Fertilizer $</span>2000&quot; = &quot;red&quot;, 
                                 &quot;Fertilizer $3000&quot; = &quot;green&quot;)) +
    labs(title = &quot;Crop Yield Forecast - Fertilizer Cost Scenarios&quot;,
         x = &quot;Date&quot;,
         y = &quot;Crop Yield&quot;,
         color = &quot;Scenario&quot;,
         fill = &quot;Forecast Interval&quot;) +
    theme_minimal() +
    theme(legend.position = &quot;bottom&quot;)
</code></pre>
<p><a href="https://i.sstatic.net/CbJHdZqr.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/CbJHdZqr.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2024-08-16T06:24:44.830Z
I'd have to read it more carefully and slowly to have any chance of saying anything useful. But just two quick things ( If I have time, I'd like to come back to your question because nerlove is a hero of mine ).  1) Check out the literature by Marc Nerlove. He's done a lot of working in econometrics and agriculture. Check out his papers because he may hae done something atleast similar to what you are doing.  If you are not familar with him, he's kind of a legend in time-series econometrics and specifically as it applies to agriculture.
mark leeds
CC BY-SA 4.0
2024-08-16T06:30:20.007Z
My other comment is to try one lag of each variable rather than two. Using a second lag of the same variable, for me, doesn't really make that much sense from an impulse response viewpoint. I could explain more why I feel this way but, for now, it's probably just best to try it and see if it makes much of a difference. I try to apply the principal of parsimony  when building these types of models.  Last thing:  I'm not familiar with auto.arima in R but you havw an ARDL so there are no moving average terms.  So, are you sure that auto.arima fits that ARDL ?
mark leeds
CC BY-SA 4.0
2024-08-16T06:41:37.743Z
Sadly, Dr. Nerlove died quite recently.  In case that you are interested, his CV contains a listing of all of his papers so it might be useful to check out that list for relevance.  https://www.iza.org/people/fellows/624/marc-l-nerlove
mark leeds
CC BY-SA 4.0
2024-08-16T14:22:04.563Z
@ Mark Leeds: Thank you so much for your reply! My condolences to Dr. Nerlove, may his family find courage and strength during this difficult time.
stats_noob
CC BY-SA 4.0
2024-08-16T14:32:14.013Z
Sure, I can look at one lag. Can you please explain why you think " a second lang doesnt really make much sense from an impulse response viewpoint"?
stats_noob
CC BY-SA 4.0
2024-08-16T17:03:40.183Z
Thanks for your kind words. I only knew Dr. Nerlove through email but he was quite a generous person. I'm writing a short answer to your question below.
mark leeds
<econometrics>
First Time Working on a Economics Problem
stats_noob
# First Time Working on a Economics Problem

(Asked by: username_0 on 2024-08-16)

I am working on this economics problem and feel a bit lost. I was hoping someone could help me out.

* Suppose there is a farm in which there are two variables: Crop Yield and Fertilizer Cost
* I am interesting in making a time series model to model the relationship between Crop Yield and Fertilizer Cost (I.e how much Fertilizer I choose to buy, assuming I fix the amount of money I am willing to spend on Fertilizer as constant)
* Specifically, I believe that Crop Yield is influenced by the (lag-1) and (lag-2) values of Crop Yield and Fertilizer Cost.
* After doing some research, it seems to me that an ARIMAX model is a good choice for this problem.
* In the end, I want to answer the following questions:
  + **What is the (estimated) impact of each extra dollar of fertilizer cost on crop yield?**
  + **Suppose I have 3 options over the next 12 months : Invest 1000 dollars each month on fertilizer, Invest 2000 dollars each month on fertilizer, Invest 3000$ each month on fertilizer. Which of these options will result in the largest (estimated) cumulative crop yield over 12 months?**

I tried to express this in mathematical terms:

Let $ Y\_t $ represent the crop yield at time $ t $, and let $ F\_t $ represent the fertilizer cost at time $ t $. We assume that the price of fertilizer remains constant, so variations in $ F\_t $ reflect the quantity of fertilizer purchased.

The crop yield $ Y\_t $ is influenced by:

* The crop yield at lag-1 ($ Y\_{t-1} $) and lag-2 ($ Y\_{t-2} $),
* The fertilizer cost at lag-1 ($ F\_{t-1} $) and lag-2 ($ F\_{t-2} $).

The relationship between crop yield and fertilizer cost can be modeled using an ARIMAX model as follows:

$$
Y\_t = \alpha\_0 + \alpha\_1 Y\_{t-1} + \alpha\_2 Y\_{t-2} + \beta\_1 F\_{t-1} + \beta\_2 F\_{t-2} + \epsilon\_t
$$

where: $ \alpha\_0 $ is the intercept term. $ \alpha\_1, \alpha\_2 $ are the coefficients for the lagged crop yields, $ \beta\_1, \beta\_2 $ are the coefficients for the lagged fertilizer costs, $ \epsilon\_t $ is the error term.

**Thus, it seems to me that $\beta$ would represent the estimated impact of fertilizer on crop yield ... and I could add the estimated 12 month forecasts for each option to find out the net benefit? Is my understanding of Economics correct?**

Can someone please tell me if I am on the right track?

Thanks!

* Note: I tried to make an R simulation to represent my problem.

First I simulate data:

```
library(forecast)
library(ggplot2)
library(dplyr)
library(tsibble)

set.seed(123)
n <- 100
dates <- seq(as.Date("2020-01-01"), by = "month", length.out = n)
fertilizer <- runif(n, 500, 1500)
seasonal_effect <- sin(2 * pi * (1:n) / 12) * 200
trend <- seq(0, 500, length.out = n)
crop_yield <- 0.5 * fertilizer + seasonal_effect + trend + arima.sim(list(ar = c(0.7, 0.2)), n) + rnorm(n, 0, 50)

data <- tibble(
    date = dates,
    crop_yield = crop_yield,
    fertilizer = fertilizer
) %>%
    as_tsibble(index = date)

```

Then, I fit the model:

```
model <- auto.arima(data$crop_yield, xreg = cbind(data$fertilizer))

```

From there, I prepared the data for the forecasting:

```
future_dates <- seq(from = max(dates), by = "month", length.out = 13)

generate_forecast <- function(cost) {
    future_fertilizer <- c(tail(data$fertilizer, 1), rep(cost, 12))
    forecast_result <- forecast(model, xreg = future_fertilizer, h = 13)
    
    list(
        date = future_dates,
        crop_yield = as.numeric(forecast_result$mean),
        fertilizer = future_fertilizer,
        lower = as.numeric(forecast_result$lower[, 2]),
        upper = as.numeric(forecast_result$upper[, 2]),
        scenario = paste0("Fertilizer $", cost)
    )
}

```

Then, I fit the forecasts and plotted them:

```
forecast_1000 <- generate_forecast(1000)
forecast_2000 <- generate_forecast(2000)
forecast_3000 <- generate_forecast(3000)

ggplot() +
    geom_line(data = data, aes(x = date, y = crop_yield, color = "Historical")) +
    geom_point(data = data, aes(x = date, y = crop_yield, color = "Historical")) +
    geom_line(data = tibble(date = forecast_1000$date, crop_yield = forecast_1000$crop_yield), aes(x = date, y = crop_yield, color = "Fertilizer $1000")) +
    geom_line(data = tibble(date = forecast_2000$date, crop_yield = forecast_2000$crop_yield), aes(x = date, y = crop_yield, color = "Fertilizer $2000")) +
    geom_line(data = tibble(date = forecast_3000$date, crop_yield = forecast_3000$crop_yield), aes(x = date, y = crop_yield, color = "Fertilizer $3000")) +
    geom_ribbon(data = tibble(date = forecast_1000$date, ymin = forecast_1000$lower, ymax = forecast_1000$upper), aes(x = date, ymin = ymin, ymax = ymax, fill = "Fertilizer $1000"), alpha = 0.2) +
    geom_ribbon(data = tibble(date = forecast_2000$date, ymin = forecast_2000$lower, ymax = forecast_2000$upper), aes(x = date, ymin = ymin, ymax = ymax, fill = "Fertilizer $2000"), alpha = 0.2) +
    geom_ribbon(data = tibble(date = forecast_3000$date, ymin = forecast_3000$lower, ymax = forecast_3000$upper), aes(x = date, ymin = ymin, ymax = ymax, fill = "Fertilizer $3000"), alpha = 0.2) +
    scale_color_manual(values = c("Historical" = "black", 
                                  "Fertilizer $1000" = "blue", 
                                  "Fertilizer $2000" = "red", 
                                  "Fertilizer $3000" = "green")) +
    scale_fill_manual(values = c("Fertilizer $1000" = "blue", 
                                 "Fertilizer $2000" = "red", 
                                 "Fertilizer $3000" = "green")) +
    labs(title = "Crop Yield Forecast - Fertilizer Cost Scenarios",
         x = "Date",
         y = "Crop Yield",
         color = "Scenario",
         fill = "Forecast Interval") +
    theme_minimal() +
    theme(legend.position = "bottom")

```

[![enter image description here](https://i.sstatic.net/CbJHdZqr.png)](https://i.sstatic.net/CbJHdZqr.png)

**Tags:** econometrics

> **Comments:**
> **username_1**: I'd have to read it more carefully and slowly to have any chance of saying anything useful. But just two quick things ( If I have time, I'd like to come back to your question because nerlove is a hero of mine ).  1) Check out the literature by Marc Nerlove. He's done a lot of working in econometrics and agriculture. Check out his papers because he may hae done something atleast similar to what you are doing.  If you are not familar with him, he's kind of a legend in time-series econometrics and specifically as it applies to agriculture.  
> **username_1**: My other comment is to try one lag of each variable rather than two. Using a second lag of the same variable, for me, doesn't really make that much sense from an impulse response viewpoint. I could explain more why I feel this way but, for now, it's probably just best to try it and see if it makes much of a difference. I try to apply the principal of parsimony  when building these types of models.  Last thing:  I'm not familiar with auto.arima in R but you havw an ARDL so there are no moving average terms.  So, are you sure that auto.arima fits that ARDL ?  
> **username_1**: Sadly, Dr. Nerlove died quite recently.  In case that you are interested, his CV contains a listing of all of his papers so it might be useful to check out that list for relevance.  https://www.iza.org/people/fellows/624/marc-l-nerlove  
> **username_0**: @ Mark username_1: Thank you so much for your reply! My condolences to Dr. Nerlove, may his family find courage and strength during this difficult time.  
> **username_0**: Sure, I can look at one lag. Can you please explain why you think " a second lang doesnt really make much sense from an impulse response viewpoint"?  
> **username_1**: Thanks for your kind words. I only knew Dr. Nerlove through email but he was quite a generous person. I'm writing a short answer to your question below.  

---

## 2 answers

### Answer by username_1 on 2024-08-19. Score: 2 | Accepted answer

Hi stats\_noob: I read through your question more carefully so I'm putting some questions in another answer. It's definitely NOT answer to your problem but contains some comments and confusions.

1a) I like the idea of using the predictions ( you call them expected values and they can be thought of this way. I like to think of them as predictions ) to make your decisions. I don't think there are constraints in your problem but I'm not totally sure. If there are, then this type of solution may not be feasible and a dynamic optimization might be a better way to go. This leads to 1b).

1b) Dynamic optimization is used when one has constraints in their problem along with some kind of model. I hope you don't have constraints because it's kind of a macro-econ field in itself and somewhat involved. I can provide some intro material on it if there are constraints. This leads into 2).

2 ) At first I thought that the amount of money that you could spend on fertilizer was constrained and that you would need to spend more when fertilizer cost was high and less when fertilizer cost was low ?
But, after reading the complete question, I don't think that's the case.

But then I asked: "Why is the fertilizer cost changing ?". I'm not clear on that which kind of makes it difficult to provide any insights. Also, how does fertilizer cost act as a predictor for crop yield ? If the fertilizer cost changes, should one expect that to increase or decrease the crop yield ? what are the units of crop yield and fertilizer cost ? Maybe fertilizer cost is not the $F\_t$ in the model. You used the term "fertilizer" so maybe $F\_t$ denotes the amount of ferilizer used at time $t$.

3 ) Also, I'm sure that auto.arima is quite a useful R function for estimation but, for your particular model, it's easier to just make a data.frame out of the 4 predictors and then use that in a call to lm. You're clearly very good with R so I think you can do that without my help. If you want some, let me know. Once you do that, then testing various models is just a straightforward call to lm and then you don't have to worry about all of these arima coefficients lying around. Personally, I still don't know what's under the hood of auto.arima but I know with certainty that you don't need it.

4 ) I realize that your R code is just a template-simulation ( not the real data etc. It's beautiful R code by the way. I've never used tidyverse or ggplot but I can still tell that you know your stuff ) I'm not clear on why you're predicting 13 steps ahead ? Shouldn't you be predicting 1 step ahead, constructing a forecast and then re-estimating the model for the next period so that you can make another forecast ? I might be confused about this also ?

So, if you can address my issues-questions, maybe I can give some useful feedback instead of confusion !!!!!!!

P.S. : 1 R tip: Don't call your data, "data". That can be dangerous because there was a function in Splus called data().
I think there's one in R also so doing that can cause confusion or even more serious problems than that.

#===============================================================
ADDENDUM: ADDING RCODE THAT SHOWS HOW TO USE lm TO ESTIMATE THE MODEL DESCRIBED IN THE ORIGINAL QUESTION. DONT KNOW WHY I HAVE TO SKIP LINES. IF I DO NOT, THINGS GET UGLY.

```
library(zoo)

set.seed(123)

n <- 10

dates <- seq(as.Date("2020-01-01"), by = "month", length.out = n)

fertcost <- runif(n, 1000, 1500)

cropyield <- runif(n, 1000, 1500)

df <- data.frame(dates, fertcost, cropyield)

mydata <- read.zoo(df)

mylaggeddata <- merge(Ytee = mydata[,2], Yteem1 = lag(mydata$cropyield, -1,), Yteem2= lag(mydata$cropyield, -2), Fteem1= lag(mydata$fertcost, -1,),
                      Fteem2 = lag(mydata$fertcost, -2))

mylaggeddata <- as.data.frame(mylaggeddata)

lm(data = mylaggeddata, as.formula(Ytee ~ Yteem1 + Yteem2 + Fteem1 + Fteem2))

```

ML: 08-20-2024
#====================================================================

Hi stats\_noob: I hope that the R code helped. Notice that, once the observations are aligned correctly, your model is no longer a time series regression model in the sense that time has been eliminated. Therefore, you just have a regression model which can be estimated using lm. Of course, if the error terms was not white noise and there was a more complicated error term structure, then you couldn't use lm.

I read through your answers to my question and I like the idea of forecasting each month out for 12 months and adding them up. The only problem with that approach is that you need the values of $F\_{t-1}$ and $F\_{t-2}$ which you won't have as you go forward. So, that approach is probably not feasible ?

The other confusion I had was that I'm still not sure if $F\_t$ is amount of fertilizer purchased or fertilizer cost ? In either case, you need to justify their ability as predictors of yield because I'm not sure of their predictive ability ? Fertilizer cost would make sense if you were spending more for a "more powerful" fertilizer but, other than that, I don't see their predictive ability ?

I hope the R code and my questions help ?

> **Comments:**
> **username_0**: @ Mark username_1: Thank you so much for all your efforts and kind words!  
> **username_0**: About 4) . I am predicting 13 steps ahead because I think this is a realistic way to test models. Farmers usually need to forecast 12 months ahead and allow 12 months worth of uncertainty into their decision making. Of course predict 1 month, re-estimate, predict 1 month, re-estimate would give you much better results... but this would not be realistic in my opinion. In my opinion, a farmer might need to buy 12 months worth of fertilizer at once ... regardless of whether half way through the year the farmer's decision turned out to be sub-optimal  
> **username_0**: it reminds me of a situation in college. one time, a chemistry professor said there is a midterm and a final. exam if you skip the midterm, the final exam will be worth the entire course mark. I innocently told the professor : "on the day of the midterm, can I sit inside, quickly flip through the midterm and see how confident I am .. and then decide if I want to write it or leave"? Incredulous, the prof replied "this is a chemistry midterm, not a buffet! where do you students come up with these ideas? if you spent half the time studying as you do coming up with these ideas you would get an A+"  
> **username_0**: The same way, I imagine that the farmer cant "take a peak at the farm conditions each month and re-update his fertilizer purchases ... return excess amounts and place more fertilizer orders if necessary". The farmer is forced to make a decision for 12 months at once  
> **username_0**: About 2) I wanted to make it that in the past, the farmer had spent varying amounts on fertilizer each month ... but when planning for the future, the farmer decides to split the year's cost over 12 months  
> **username_0**: I would be very curious to hear your lm() approach in R. If you have time, can you please write it?  
> **username_0**: BTW - have you heard of linear state space models?  
> **username_1**: let me print out your comments and get back to you in a day or two. I wrote something in R ( I read your things quickly ) where I do an lm on artificial data. The sad thing is that I don't know how to put it in that window that you put yours in. Can you explain how to do that ? The code is short and simple. You'll see what's going on right away. Thanks.  
> **username_1**: Yes, I have some experience with state space models. But I always try to figure out which variables might be helpful and possible relationships between them BEFORE I decide on a methodology. That's kind of secondary, atleast from  my viewpoint. Others may feel differently. It's all good. But yes, a kalman filter is useful when you don't necessarily believe that your predictors are constant. Basically, the predictors vary over time, possibly. Other than that reason, I see no use for it with your current model.  
> **username_1**: If the window creation thing is complicated, I can just put the R code in another answer ? But that would be my third answer to the same question !!!!! I imagine that's discouraged by whomever are the kind people who run this list ? Shout out to them. It's a wonderful list with many generous and kind people.  
> **username_1**: I managed to make the R code to use lm an addendum in the second answer. Only negative is that I had to take out my comments to get it to look decent. If any of the R code is not clear, just let me know. I have comments in the actual code.  
> **username_1**: @username_0: I read through your answers and put another addendum below the addendum containing the R-code. I hope it helps.  
> **username_2**: I have edited the answer. Use curly-brackets symbol '{}' for R-code.  
> **username_1**: Thanks username_2: I'll go over it to see what you did. It's much appreciated.  
> **username_1**: username_2: I went into edit mode to see what you did but I can't find the {}. I must have woken up too early this morning. It definitely looks much better so thanks but I can't figure out what you did to do that. uuuugh.  

---

### Answer by username_1 on 2024-08-16. Score: 2

Let me just give a short example as to why two lags of $y\_t$ might be "redundant" if you will.

Suppose you just had an AR(1), so

$$y\_t = \beta y\_{t-1} + \epsilon\_t$$

Then, what is this model REALLY saying. Converting it using the lag operator, one obtains:

(1)
$$y\_t = \sum\_{i=0}^{\infty} \beta^{i} \epsilon\_{t-i} $$

So, assume a shock of $\epsilon\_t$ at time $t$ ( this is referred to as the impulse in the econometrics literature ) and assume no shocks after that. Then this model implies that the impulse response of a shock is an exponentially smoothed version of that same shock. So, let us call that shock $\epsilon\_t$ at time $t$ as $\epsilon^{\*}$ and assume that there are no shocks after that.

Then, based on (1), first, at time $t$ the effect of the $\epsilon^{\*}$ shock is just $\epsilon^{\*}$. Then, at time $t+1$, the shock's effect get decreased to $\beta \epsilon^{\*}$. Then at time $t+2$, this same shock has effect of $\beta^2 \epsilon^{\*}$ and so on and so forth until the effect of the shock totally dissipates and eventually becomes zero.

Now, suppose the model had been an AR(2), rather than an AR(1). So, in this case, one has:

$$y\_t = \beta y\_{t-1} + \theta y\_{t-2} + \epsilon\_t $$

So, going through a similar derivation one obtains $$y\_t = \sum\_{i=0}^{\infty}\theta^{i} \epsilon\_{t-i-1}$$

So not only do we have this one set of "effects" due to the
shock $\epsilon^{\*}$ at time $t$. We then have this second set of exponentially smoothed effects due to the same $\epsilon^{\*}$. This set starts at time $t+1$ and result in sequential effects that are exponentially smoothed over time. The difference between this second set of shocks and the first set is that they start at time $t+1$ and have a different smoothing coefficient with the coefficient being $\theta$ rather than $\beta$.

My point about lagging once rather than twice is that I find this interpretation for the second set of effects difficult to justify when the shock they originate from is the same shock just lagged by one more time unit. Why would there be this "similar" set of effects except lagged by one time unit ? The first set already takes care of the time lagging effect decreasing over time so there probably shouldn't be the need for a second set ? It seems like too much overlapping to me unless there was some known reason for it.

I provided the argument for the AR(1) but a similar argument could be used for the ARDL that you're using. The $F\_t$ are assumed exogenous so one cannot argue using the $\epsilon\_t$s like we did in the AR examples above but still a similar argument holds: The effects of the $F\_t$ coefficients also dissipate over time due to the auto-regressive-ness of the model so their set of effects are also exponentially smoothed over time. This is why I would argue against a second lag for the $F\_t$ also.

When I was reading about ARDLs ( I think they are synonymous with what is referred to by statisticians as ARIMAX but I'm not certain ), I found that there is a lot of good political science literature on distributed lag models. This was somewhat surprising and it helped me a lot. The authors that I highly recommend are Luke Keele, Nathaniel Beck and Suzanne Deboef. If you can't find anything connected to them and distributed lags, let me know and I can send some links. Some of their papers were real eye-openers for me and even better than most of explanations in the econometrics texts. All the best.

ADDENDUM TO ORIGINAL POST
#================================================================

Notice that the exponentially smoothed "set" of effects arise
ONLY WHEN THERE IS A LAGGED AUTO-REGRESSIVE term, in this case,
$y\_{t-1}$ and $y\_{t-2}$. If there were no lagged auto-regressive
terms then, any shock or exogenous effect could be viewed as an impulse but it would be a one time impulse at time $t$ and then zero after that. So, this "lagged effect" concept from the impulse response would not exist without the existence of atleast one of the lagged auto-regressive terms, $y\_{t-1}$ or $y\_{t-2}$.

SECOND ADDENDUM TO ORIGINAL POST
#==========================================

Note that when explaining the two "sets of effects" caused by having two lagged versions of the same variable, I simplified things are far as the two sets being seperate and additive. The effect from the first lag is being fed into the second lag so the effects are definitely not simply seperate and additive. One would have use the polynomial equivalence approach used in transfer function modelling to figure out the exact sequence of effects due to the impulse response. Nevertheless, even though the way of thinking about the two sets being "additive" is not correct, the argument still holds: "Why would an impulse response of $\epsilon^{\*}$ cause such a complex relationship for the effects of that response over time. It's possible but not likely.

> **Comments:**
> **username_0**: @ Mark username_1: thank you so much for your detailed answer! Just a question - does the overall idea sound right? we can predict the expected values - and then based on the predicted expected values conduct the cost benefit analysis?  
> **username_1**: Hi: Youre welcome. But that last part is really important as far as understanding that you can just take two seperate effects and look at them individually ( like I did ) because what you really have are two filters in series so  SYSTEM RESPONSE would need to be figured out.  So I simplified things terribly when I made a feeble attempt to explain what's being assumed about the impulse response when an AR(2) is used.  
> **username_1**: Hi username_0: Let me read it again over the next couple of days because I really didn't read it the way it should be read. I just focused on the estimation part. So, let me get back to you. It sounds like an important question so I will just say whatever I would do if I was doing what you're doing. Others may have different ways ?  I've realized in the last 10 or so years that the amount of overlap in different fields is enormous but even the overlap varies because it depends on what is emphasized. That's the most interesting thing: There's variation in the overlap itself !!!!!!  
> **username_1**: first sentence of next to last comment is "you can't just take". To figure out the overall response given the AR(2) you have to compute the transfer function. This sort of thing falls more under the topic of "system identification". It is also covered in "signals and systems" type classes. For whatever reason, it's  more emphasized in digital signal processing. But we're getting off track so I'll leave it at that.  
> **username_1**: username_0: here's a link to Hyndman's text where he explains how the auto.arima model works. I'm just sending it now so that I don't forget. It would be useful to find an example where auto.arima is used to fit a model like yours rather than a reg model with arima errors. He likes the latter approach which has equivalences to the  ARDL models which are more in line with your notation-setup.  I haven't done this sort of thing ( fit ARDL models ) in a while but the systemfit package might be a better option for estimation of your model.  https://otexts.com/fpp2/dynamic.html  

---

