System:

Y_t = C_{t} + I_t + G_t

C_t = a + \alpha Y_{t-1}

I_t = \beta (Y_{t} - Y_{t-1})

G_t = 0\times 1\{t\leq 14\} + 100\times 1\{t\geq 15\}

Parameters: a=50, \alpha =.75, \beta =4

Initial Condition: Y_{0}

System of 4 equations in 4 variables Y_{t}, C_{t}, I_{t}, G_{t}

A solution gives each variable as a function of `t`

.

G_{t} already given (exogenously).

**Solution 1 (Closed form)**

Plug the 2nd/3rd equation into the 1st:

Y_t = (a + \alpha Y_{t-1}) + (\beta Y_{t} - \beta Y_{t-1}) + G_t

\Leftrightarrow (1-\beta)Y_t = a + (\alpha -\beta) Y_{t-1} + G_t

\Leftrightarrow Y_t = \frac{a}{(1-\beta)} + \frac{(\alpha -\beta)}{(1-\beta)} Y_{t-1} + \frac{1}{(1-\beta)}G_t , Y_0 given is a 1st order difference equation, solve for Y_{t}.

Next, plug Y_{t} into the 2nd/3rd equations to get C_{t}, I_{t}

Hopefully `Symbolics.jl`

will one day solve systems of difference equations @ChrisRackauckas.

For now Mathematica:

**Solution 2 (numerical in Julia)**

```
using Plots;
#
a=50.0; α=.75; β=4.0; Y0=0.0;
G(t)= (t <= 14) ? 0.0 : 100.0
Y_rec(Ym, t; α=α, β=β, a=a) = (a/(1-β)) + ((α-β)/(1-β))*Ym + (1/(1-β))*G(t)
I_rec(Y, Ym; β=β) = β*(Y-Ym)
C_rec(Ym; a=a, α=α) = a + α*Ym
# Simulate
Δt = 1.0; T = 250; time = 0.0:Δt:T
G_sim, Y_sim, I_sim, C_sim = [zeros(length(time),1) for i in 1:4]
Y_sim[1] = Y0
#
for i in 1:(length(time)-1)
G_sim[i+1] = G(i)
Y_sim[i+1] = Y_rec(Y_sim[i], i)
I_sim[i+1] = I_rec(Y_sim[i+1], Y_sim[i])
C_sim[i+1] = C_rec(Y_sim[i])
end
#
ix = 2:21 #(100) # time=0, to time=23
p1 = plot(size=(400,400),legend=:bottomleft)
plot!(time[ix], G_sim[ix], lab = "G")
plot!(time[ix], Y_sim[ix], lab = "Y")
plot!(time[ix], I_sim[ix], lab = "I")
plot!(time[ix], C_sim[ix], lab = "C")
```

If we change one parameter \beta =0.4 instead of \beta=4, things look more as expected:

As we lengthen the time horizon, variables approach Steady State:

@josec Based on the simulation, I believe there was a typo somewhere in your description of the problem (or parameter value, or initial value). Can you please post the original source?