This page gives an example of code you can use to numerically solve the Solow model under various conditions and use that to plot figures. It's a Quarto document that you can edit in RStudio. You can [download the original Quarto doc here](quarto-solow.qmd), and once you open it in RStudio you can "Render" the document which will produce a webpage (in html) like this one. Or if all you want to do is copy the actual R *code* that is being run to create the figures, you can just cut/paste the code itself. 

## Solving the Solow
This first block of code is a function that solves the Solow model. The function takes parameters like the capital accumulation rate, sI, or the population growth rate, gL. It also takes initial conditions for capital, labor, and productivity. The last, and maybe most important, parameter is the time periods to solve for.

You can use the code "as-is". If you want to understand what's going on in the actual calculation, look at the Study Guide for the "An exact solution to the Solow model". This key code is the line starting "ky_time <-". Solving the Solow is basically about solving for the exact K/Y ratio at each point in time given the parameters.

```{r}
#| code-fold: false

# You need this plotly library to do the figures. This just absorbs any messages about loading it when you render.
suppressMessages(library(plotly))

f.solow <- function(alpha = .3, si=.2, delta=.05, ga=.02, gl=.01, k0=1, a0=1, l0=1, t=1) {
  kal_init <- k0/(a0*l0) # initial K/AL ratio
  ky_init <- kal_init^(1-alpha) # initial K/Y ratio
  ky_ss <- si/(delta+ga+gl) # K/Y steady state
  kal_ss <- ky_ss^(1/(1-alpha)) # K/AL steady state
  gK_init <- si/ky_init - delta # initial gK
  gy_init <- (1-alpha)*(gK_init - ga - gl) + ga # initial gy
  
  # Manage time periods passed
  u <- tryCatch(
    {u = t - min(t)}, # reset any list of time periods to start at zero, for convergence calcs
    error=function(cond) { # if above fails, just take t as-is (usually a single number)
      return(t)
    }
  )
  
  conv_term <- (1-alpha)*(delta+ga+gl) # convergence term for exact solution
  ky_time <- ky_ss*(1-exp(-1*conv_term*u)) + ky_init*exp(-1*conv_term*u) # K/Y at each u
  kal_time <- ky_time^(1/(1-alpha)) # K/AL at each t
  lny_time <- (alpha/(1-alpha))*log(ky_time) + log(a0) + ga*t # Log GDP pc at each t
  lny_bgp <- (alpha/(1-alpha))*log(ky_ss) + log(a0) + ga*t # BGP of GDP pc at each t
  gK_time <- si/ky_time - delta # growth rate of K at each t
  gy_time <- (1-alpha)*(gK_time - ga - gl) + ga # growth rate of GDP pc at each t 
 
  # create list containing all results
  results <- list("time" = t, "KAL" = kal_time, "KY" = ky_time, "lny" = lny_time, 
                 "lny_bgp" = lny_bgp, "gK" = gK_time, "gy" = gy_time, "KY_ss" = ky_ss,
                 "KAL_ss" = kal_ss, "KY_init" = ky_init, "KAL_init" = kal_init, 
                 "gK_init" = gK_init, "gy_init" = gy_init,
                 "alpha" = alpha, "si" = si, "delta" = delta, "gA" = ga, "gL" = gl)
  solved <- as.data.frame(results) # convert to data frame to pass back

  return(solved) # return all results
}
```

## Create basic figures
Let's use this function to map out some scenarios. This starts by telling R that we want to plot the function from year 1 to year 60. Then it calls the Solow function to calculate all the outcomes over those years. An important point here is that you don't have to specify everything for the function. You can just specify the thing you think is interesting, letting the function set the other parameters.


```{r}
#| code-fold: false
years = seq(from = 1,to = 60, by = 1)
s= f.solow(t = years)

```

Given that solution we can plot figures of what is going on

```{r}
#| label: fig-level
#| fig-cap: "Log GDP per capita over time"
#| warning: false

# Plot figure of scenarios and their BGP's
fig1 <- plot_ly(s, x = ~time, y = ~lny, type = 'scatter', mode = 'lines', name='Log GDP pc')
fig1 <- add_trace(fig1, x = ~time, y = ~lny_bgp, type = 'scatter', mode = 'lines', name='BGP')

fig1
```

Then do a figure of the growth rate of GDP per capita. Note how this matches up to the path for log GDP per capita.

```{r}
#| label: fig-growth
#| fig-cap: "Growth rate over time"
#| warning: false

fig2 <- plot_ly(s, x = ~time, y = ~gy, type = 'scatter', mode = 'lines', name='Growth rate')
fig2 <- add_trace(fig2, x = ~time, y = ~gA, type = 'scatter', mode = 'lines', name='gA')
fig2 <- layout(fig2, yaxis = list(title="Growth rate",range=c(0,.1)))

fig2
```

Last, you can visualize the path of K/AL and the growth rate of capital during this transition. 

```{r}
#| label: fig-kal
#| fig-cap: "K/AL dynamics"
#| warning: false

fig3 <- plot_ly(s, x = ~KAL, y = ~gK, type = 'scatter', mode = 'lines+markers', name='Capital growth rate')
fig3 <- add_segments(fig3, data = s, x = 0, xend = 4, y = ~(gA+gL), yend = ~(gA+gL), name='gA + gL', line = list(dash = "dash"))
fig3 <- add_segments(fig3, data = s, x = ~KAL_ss, xend = ~KAL_ss, y = 0, yend = .15, name = 'K/AL steady state', line = list(dash = "dash"))
fig3 <- add_segments(fig3, data = s, x = ~KAL_init, xend = ~KAL_init, y = 0, yend = .15, name = 'Initial K/AL', line = list(dash = "dash"))
fig3 <- layout(fig3, yaxis = list(title="Capital growth rate",range=c(0,.15)))
fig3 <- layout(fig3, xaxis = list(title="K/AL ratio",range=c(0,4)))

fig3

```

## Show a shift in the capital accumulation rate
We can do something different. Here's an example of starting the model *at* steady state, and then running the code with a relatively high savings rate to see how the economy shifts to a higher BGP.

```{r}
#| code-fold: false
years = seq(from = 1,to = 60, by = 1) # run for this many periods
b = f.solow(t = years) # run the model at baseline to get some parameters
```

From that model the value KAL_ss contains the SS K/Al ratio. With A = 1 and L =1 we're going to set the initial capital stock K0 to the KAL_ss value so the model starts at steady state.

```{r}
#| code-fold: false

s= f.solow(si = .4, t = years, k0 = b$KAL_ss) # solve using the higher capital accum rate
o= f.solow(t=years, k0 = b$KAL_ss)

```

Now let's go back and plot the path of GDP per capita. This figure includes a plot from the baseline of the *old* baseline BGP for comparison.

```{r}
#| label: fig-compare-bgp
#| fig-cap: "Log GDP per capita over time"
#| code-fold: false
fig1 <- plot_ly(s, x = ~time, y = ~lny, type = 'scatter', mode = 'lines', name='Log GDP pc')
fig1 <- add_trace(fig1, x = ~time, y = ~lny_bgp, type = 'scatter', mode = 'lines', name='BGP')
fig1 <- add_trace(fig1, data = o, x=~time, y = ~lny_bgp, type = 'scatter', mode = 'lines', name='Old BGP' )

fig1
```

We can do something similar for the plot of growth rates over time.

```{r}
#| label: fig-compare-rate
#| fig-cap: "Growth rate over time"
#| warning: false

fig2 <- plot_ly(s, x = ~time, y = ~gy, type = 'scatter', mode = 'lines', name='Growth rate')
fig2 <- add_trace(fig2, data = o, x = ~time, y = ~gy, type = 'scatter', mode = 'lines', name='Baseline')
fig2 <- layout(fig2, yaxis = list(title="Growth rate",range=c(0,.1)))

fig2
```

To try and plot the K/AL dynamics in the baseline and then in the alternative scenario, we want to run another Solow model to gather info on the growth rate of K under the baseline at various K/AL ratios.

```{r}
#| code-fold: false

kal = seq(from = 0,to = 10, by = .1) # sequence of values for K
p= f.solow(k0 = kal) # solve model for all values of K, which is =K/AL given A=1, L=1

```

Now we can plot two relationships of gK to K/AL. The figure shows that the whole curve of capital growth rates, gK, shifted up when sI shifted to 0.4 (the blue line is above the orange). 

```{r}
#| label: fig-kal-comp
#| fig-cap: "K/AL dynamics versus baseline"
#| warning: false

fig3 <- plot_ly(s, x = ~KAL, y = ~gK, type = 'scatter', mode = 'lines+markers', name='New capital growth')
fig3 <- add_trace(fig3, data = p, x = ~KAL, y = ~gK, type = 'scatter', mode = 'lines+markers', name = 'Baseline capital growth')
fig3 <- add_segments(fig3, data = s, x = 0, xend = 12, y = ~(gA+gL), yend = ~(gA+gL), name='gA + gL', line = list(dash = "dash"))
fig3 <- add_segments(fig3, data = s, x = ~KAL_ss, xend = ~KAL_ss, y = 0, yend = .15, name = 'New steady state', line = list(dash = "dash"))
fig3 <- add_segments(fig3, data = p, x = ~KAL_ss, xend = ~KAL_ss, y = 0, yend = .15, name = 'Baseline steady state', line = list(dash = "dash"))
fig3 <- layout(fig3, yaxis = list(title="Capital growth rate",range=c(-.04,.15)))
fig3 <- layout(fig3, xaxis = list(title="K/AL ratio",range=c(0,12)))

fig3

```


