Replicating Solow Growth Model in R

This page will illustrate the replication process from the example of the Solow Growth Model from the book “An introduction to R for quantitative economics. SpringerBriefs in Economics.” of Vikram Dayal.

This is done entirely for educational purposes. All credits to the idea, source code belong to Vikram Dayal. The citation for the book is:

Dayal, V. (2015). An introduction to R for quantitative economics. SpringerBriefs in Economics.https://link.springer.com/book/10.1007/978-81-322-2340-5

The Author of the post: John Riveros

In this Markdown we’re going to mimic the programming of the Solow model of Vikram Dayal

Libraries to use

We will use the mosaic package to create the functions and plots. If it is not installed before, uncomment the next command and then load up the library.

#install.packages("mosaicCalc")
library(mosaicCalc)
 #install.packages("mosaic")
library(mosaic)
#install.packages("latticeExtra")
library(latticeExtra)

Production function

We start by defining the production function as

\[ Y = A K^a L^{1-a} \] Here \(Y\) refers the output of the economy, \(A\) is a techonological parameter, \(K\) is the aggregate capital stock, \(L\) is the aggregate labor (in time) of the economy. and \(a\) is the elasticity (or contribution) of the capital to the output \(Y\). We will not use the classic notation of \(\alpha\) as in the code, it can be more readible with just \(a\).

Parameters

It is important to note that for simplification, we can define some parameters and variables in our model. We will impose that

\[ a= 0.5 \\ A=2 \\ L=200\]

Thus the code will be based on this parameters as described next:

Y <- makeFun(A * (K^a) * (L^(1 - a)) ~ K,  #Production function
             a = 0.5,                      #elasticity
             A = 2,                        # Techonology
             L = 200)                      # Labor size

We now execute this so further we can plot the behaviour of the output of the economy when capital increases.

plotFun(Y, xlim = range(0,4000), xlab="Capital (K)", ylab="Output (Y)", 
        col = "black", lwd = 2, main = "Output - Capital Reaction")

Now we can see the dynamics, when labor is static.

Savings Curve of the economy.

We are able to define the savings of the economy as a fraction of the outcome therefore the savings function can be defined as:

\[ S = sY \] In this case, \(s\) is a marginal propensity of saving which can be interpreted as a fraction of savings from the output of the economy \(Y\).

Let’s define the savings function:

S <- makeFun(s * A * (K^a)* (L^{1-a}) ~ K, 
             s = 0.2,
             a = 0.5,                      #elasticity
             A = 2,                        # Techonology
             L = 200)                      # Labor size

Check that \(A* (K^a)* (L^{1-a})\) is equal to \(Y\). Now we can plot the savings curve.

plotFun(S, xlim = range(0,4000), xlab = "Capital (K)", ylab = "Savings (S)",
        main = "Savings curve", lwd = 2)

Depreciation.

Generally the levels of depreciation are fraction of the capital stock. We will just define the depreciation as

\[Depreciation = c K\]

We thus define the function in R as:

Depreciation <- makeFun(c * K ~ K, c = 0.1)

Then we plot the curves together

plotFun(S, xlim = range(0,4000), xlab = "Capital (K)", ylab = "S, D",
        main = "Savings and Depreciation (D)", lwd = 2)

plotFun(Depreciation, xlim = range(0,4000), xlab = "Capital (K)", 
        ylab = "Depreciation (D)", lwd = 2, col = "black", add = TRUE)

By differentiating respect to time the aggregate capital \(K\) it can be found that:

\[ dK / dt = S - Depreciation \]

Hence, the evolution of capital over time is a direct function of the aggregate levels of Savings minus de depreciation in the economy. The previous equation describes the evolution of capital (which clearly depends on the levels of aggregate income \(Y\) and the propensity to save \(s\))

We can return to our general expressions in the economy where

\[ Y = A K^a L^{1-a} \\ S = sY \\ Y = C + I \\ S = I \\ S = Y - C \\ I = dK/dt + Depreciation \\ Depreciation = cK\]

By doing some manipulation of the definition of Investment we can explore the evolution of capital over time as:

\[ dK/dt = sY - cK\] Therefore, the evolution of capital is directly a fraction \(s\) from the total income \(Y\) minus the aggregated depreciation \(cK\). Note that \(c\) is a constant fraction which is depreciated from the total aggregated capital \(K\).

By solving \(dK/dt = 0\) we can find the steady states values of aggregate capital. Hence from:

\[ 0=dK/dt = s A K^a L^{1-a} - cK \]

The \(K*\) which reflects a steady state can be found by:

Steady.state.K <- findZeros( s*A*(K^a)*(L^(1-a))-c*K ~ K,
                            s = 0.2,              #Savings rate
                            a = 0.5,              #Elasticity
                            A = 2,                #Techonology
                            L = 200,              #Labor size
                            c = 0.1               # Capital depreciation
                        )

And we can visualize it:

Steady.state.K
##      K
## 1    0
## 2 3200

Where \(K=3200\)

We can also visualize this on the previous plot:

plot(0, 0, type = "n", xlim = c(0, 4000), ylim = c(0, max(S(4000), Depreciation(4000))),
     xlab = "Capital (K)", ylab = "S, D", main = "Savings and Depreciation (D)")
curve(S(x), from = 0, to = 4000, add = TRUE, lwd = 2)
curve(Depreciation(x), from = 0, to = 4000, add = TRUE, lwd = 2, col = "darkblue")
abline(v = 3200, col = "red", lwd = 2, lty = 2)

Integration of Capital over time.

Path of aggregate capital can be calculated by solving the differential equation by integrating the previous equation. To do so. we impute the equation.

\[ dK = s*A*K^a*L^{(1-a)}, \\ K_0 = 1000 \\ \]

Here we picked up a random initial value of capital of \(K_0 = 1000\) and from here we integrate over a period of 100 years.

Solow <- integrateODE(dK ~ (s * A * (K^a) * (L^(1 - a))) -c * K ,
                      K = 1000, 
                      a = 0.5, 
                      A = 2, 
                      L = 200, s = 0.2, c = 0.1,
                      domain(t=0:100))

plotFun(Solow$K(t)~t, t.lim=range(0,100))

From this result. The dynamics of the evolution of capital can be trace down.