In this post we continue our Mankiw, Romer and Weil marathon by visualising countries’ saving and depreciation curves, derived from the Solow and Augmented Solow model, using MRW’s dataset on economic growth with the help of Plotly and Crosstalk.

Calculating alpha and beta

First, we need to estimate the model to derive empirical estimates of \(\alpha\) and \(\beta\). MRW calculate both parameters using both a restrict and un-restricted model but choose to stick with the restricted parameter estimates so we do the same.

Running the restricted regressions just involves subtracting our previous regressors from each other. Here, we use sym and !! to ensure the quasiquotation filter expects works (this is analogous to the more lengthy version we used before where we’d code each model separately):

library(purrr)
library(tidyr)

MRW_clean_restricted <- MRW_clean %>% 
  mutate(restricted_regressor = log(s) - log(n_g_d))
restricted_models_old <- c('N',
                           'I',
                           'O') %>% 
  map(sym) %>% 
  map(~filter(MRW_clean_restricted, !!. == 1)) %>% 
  map(~lm(log(`1985`) ~ restricted_regressor,
                       data = .))

Now, we collect results using collect.results defined previously and repeat the above for the Augmented model.

restricted_model_results_old <- restricted_models_old %>% 
  collect.results(., names = c("OECD Restricted Classic",
                               "Intermediate Restricted Classic",
                               "Non-Oil Restricted Classic"))

calculate.alpha <- function(estimate){
  alpha <- estimate / (1 + estimate)
  return(alpha)
}

restricted_model_results_old <- restricted_model_results_old %>% 
  mutate(alpha = ifelse((term == 'restricted_regressor'), calculate.alpha(estimate), NA))

## Augmented model
MRW_augmented_restricted <- MRW_clean_restricted %>% 
  mutate(restricted_regressor_hc = log(school / 100) - log(n_g_d))

augmented_restricted_model_old <- c("O",
                                    "I",
                                    "N") %>% 
  map(sym) %>% 
  map(~filter(MRW_augmented_restricted, !!. == 1)) %>% 
  map(~lm(log(`1985`) ~ restricted_regressor + restricted_regressor_hc,
          data = .)) %>% 
  collect.results(., names = c("OECD Restricted Augmented",
                               "Intermediate Restricted Augmented",
                               "Non-Oil Restricted Augmented"))

Calculating \(\alpha\) for the classical Solow model was a pretty simple affair as we know that the regression coefficient is equal to \(\frac{\alpha}{1-\alpha}\). The Augmented model is a little trickier since we’re solving a set of simultaneous equations. However, the equations describing \(\alpha\) and \(\beta\) are symmetric so we can just swap round y and x in the mutate function.

We use spread here because we want each estimated coefficient to have its own column - an alternative would be to use some form of apply across rows.

find.param <- function(y, x){
  param <- (x*y) / (x + x^2 + x*y)
  return(param)
}

augmented_parameters <- augmented_restricted_model_old %>%
  select(term, subset, estimate) %>% 
  filter(term != "(Intercept)") %>% 
  spread(term, estimate) %>% 
  mutate(alpha = find.param(restricted_regressor_hc, restricted_regressor),
         beta = find.param(restricted_regressor, restricted_regressor_hc),
         alpha = round(alpha, 2),
         beta = round(beta, 2)) %>% 
  select(subset, alpha, beta)

classic_parameters <- restricted_model_results_old %>% 
  select(subset, alpha) %>% 
  drop_na() %>% 
  mutate(alpha = round(alpha, 2))

estimated_parameters_old <- bind_rows(classic_parameters,
                                  augmented_parameters) %>% 
  separate(subset, c("subset",
                     "estimation_type",
                     "model_type"),
           sep = " ")
estimated_parameters_old  %>% kable(digits = 2)
subsetestimation_typemodel_typealphabeta
OECDRestrictedClassic0.60NA
IntermediateRestrictedClassic0.59NA
Non-OilRestrictedClassic0.36NA
IntermediateRestrictedAugmented0.300.29
Non-OilRestrictedAugmented0.270.31
OECDRestrictedAugmented0.370.14

One of the drawbacks of MRW’s data partitioning discussed previously is that countries aren’t partitioned into mutually exclusive groups. This is a bit of a pain because it means we effectively have multiple \(\alpha\)s and \(\beta\)s for any country which is in more than one group. For instance, the United Kingdom is included in every category by MRW so should we calculate its saving curve with an \(\alpha\) of 0.6, 0.59 or 0.36?

To overcome this we can return to the mutually exclusive partitions we created previously when we explored clustering the dataset.

Therefore, we repeat the above but with our new ‘club’ factor. In hindsight, using more functions would have been a lot quicker here:

classic_parameters_new <- c("oil",
                           "developing_no_oil",
                           "intermediate_no_oil",
                           "rich") %>%  
  map(~filter(MRW_clean_restricted, club == .)) %>% 
  map(~lm(log(`1985`) ~ restricted_regressor,
                       data = .)) %>% 
  collect.results(., names = c("Oil Restricted Classic",
                               "Developing Restricted Classic",
                               "Intermediate Restricted Classic",
                               "Rich Restricted Classic")) %>% 
  mutate(alpha = ifelse((term == 'restricted_regressor'), calculate.alpha(estimate), NA))

augmented_parameters_new <- c("oil",
                           "developing_no_oil",
                           "intermediate_no_oil",
                           "rich") %>% 
  map(~filter(MRW_augmented_restricted, club == .)) %>% 
  map(~lm(log(`1985`) ~ restricted_regressor + restricted_regressor_hc,
                       data = .)) %>% 
  collect.results(., names = c("Oil Restricted Augmented",
                               "Developing Restricted Augmented",
                               "Intermediate Restricted Augmented",
                               "Rich Restricted Augmented")) %>% 
  select(term, subset, estimate) %>% 
  filter(term != "(Intercept)") %>% 
  spread(term, estimate) %>% 
  mutate(alpha = find.param(restricted_regressor_hc, restricted_regressor),
         beta = find.param(restricted_regressor, restricted_regressor_hc),
         alpha = round(alpha, 2),
         beta = round(beta, 2)) %>% 
  select(subset, alpha, beta)

classic_parameters_new <- classic_parameters_new %>% 
  select(subset, alpha) %>% 
  drop_na() %>% 
  mutate(alpha = round(alpha, 2))

estimated_parameters_new <- bind_rows(classic_parameters_new,
                                  augmented_parameters_new) %>% 
  separate(subset, c("subset",
                     "estimation_type",
                     "model_type"),
           sep = " ")
estimated_parameters_new %>% kable(digits = 2)
subsetestimation_typemodel_typealphabeta
OilRestrictedClassic7.77NA
DevelopingRestrictedClassic0.36NA
IntermediateRestrictedClassic0.52NA
RichRestrictedClassic0.36NA
DevelopingRestrictedAugmented0.130.27
IntermediateRestrictedAugmented0.300.25
OilRestrictedAugmented0.99-1.80
RichRestrictedAugmented0.370.14

Next, we run a quick sanity check by comparing the estimates. It seems like dropping Oil is a good idea as these give pretty crazy results (somewhat predictably).

estimated_parameters_new_no_oil <- estimated_parameters_new %>% 
  filter(subset != 'Oil')
comparison_tbl <- bind_rows(estimated_parameters_new_no_oil %>%
                              select(model_type, subset, alpha, beta) %>%
                              mutate(type = 'New'),
                            estimated_parameters_old %>% select(model_type, subset, alpha, beta) %>% 
                              mutate(subset = paste0(subset, ' MRW'),
                                     type = 'Original')) 
comparison_tbl %>%
  head %>% 
  kable(digits = 2)
model_typesubsetalphabetatype
ClassicDeveloping0.36NANew
ClassicIntermediate0.52NANew
ClassicRich0.36NANew
AugmentedDeveloping0.130.27New
AugmentedIntermediate0.300.25New
AugmentedRich0.370.14New

We can visualise the differences using Plotly and plot_new_vs_old which we also defined in a previous post:

library(plotly)
comparison_plot_classic <- comparison_tbl %>% 
  filter(model_type == 'Classic') %>% 
  gather(term, estimate, alpha, beta)

plot_alpha_beta <- function(data, beta = FALSE){
  data$subset <- factor(data$subset,
                                   levels = c('Developing',
                                              'Intermediate',
                                              'Rich',
                                              'Non-Oil MRW',
                                              'Intermediate MRW',
                                              'OECD MRW'),
                                   ordered = TRUE)
  
  p <- plot_new_vs_old('alpha', data)
  if (beta == TRUE){
    q <- plot_new_vs_old('beta', data)
    p <- subplot(p, q, nrows = 2, titleY = TRUE)
  }
   return(p)
}

p <- plot_alpha_beta(comparison_plot_classic) %>% 
  layout(showlegend = FALSE,
         title = 'Calculating alpha with new and old partitions, \n Classic Solow')
p

The \(\alpha\) estimates are pretty different to MRW’s original results which is a little surprising considering the underlying regression estimates seemed to be quite similar when we plotted them in the last post. Furthermore, differentiating the \(\alpha\) formula with respect to \(x\), our regression estimates, gives us \(\frac{1}{(1+x)^2}\) which suggests that \(\alpha\) is less sensitive to any given change in \(x\) than \(x\) itself i.e. \[1 \geq\frac{1}{(1+x)^2} \quad \forall \space x\]

[Edit: This isn’t true for \(-2 < x < 0\) which invalidates the above somewhat.] This could be because the restricted regression set-up we’re employing effectively compounds all the estimated coefficient differences into one parameter.

Moving onto the Augmented model and the estimates begin to look similar to Mankiw, Romer and Weil’s:

comparison_plot_augmented <- comparison_tbl %>% 
  filter(model_type == 'Augmented') %>% 
  gather(term, estimate, alpha, beta)

p <- plot_alpha_beta(comparison_plot_augmented, beta = TRUE) %>% 
  layout(showlegend = FALSE,
         title = 'Calculating alpha and beta with new and old partitions, \n Augmented Solow')
p

Plotting classical saving and depreciation curves

Now that we have estimated our \(\alpha\) and \(\beta\) parameters, we can generate the necessary curves.

First, we create a vector of capital levels, \(k\), for each country and use left_join to include country specific variables such as s and n_g_d. Next, we use gsub to reformat the club factor entries so that we can perform another left join and add our estimated \(\alpha\)s and \(\beta\)s to the relevant countries depending on which ‘club’ they’re in.

Finally, we calculate the savings curve using the formula \(sk^\alpha\) and effective depreciation curve \((n + g + d)k\). As a (limited) test to ensure we’ve done everything correctly we calculate \(k^*\) and \(s^*\), the steady state levels of capital accumulation and corresponding savings curve point, which should hopefully correspond to where a country’s saving and depreciation curve meet.

k <- rep(seq(from = 0, to = 20, by = .5), nrow(MRW_clean))
countries <- rep(MRW_clean$country, 20/0.5 + 1)

country_tbl <- tibble('k' = k,
                      'country' = countries) %>% 
  left_join(y = MRW_clean %>% 
              select(`1985`, s, n_g_d, country, club), by = 'country')

estimated_parameters_new$club <- estimated_parameters_new$subset %>% 
  gsub('Oil', 'oil', .) %>% 
  gsub('Developing', 'developing_no_oil', .) %>% 
  gsub('Intermediate', 'intermediate_no_oil', .) %>% 
  gsub('Rich', 'rich', .)

country_tbl <- left_join(country_tbl, estimated_parameters_new %>% 
                           select(club, model_type, alpha, beta),
                         by = 'club')

find.saving.point <- function(k, alpha, s){
  point <- s*k^alpha
  return(point)
}

find.effective.depreciation <- function(k, n = 0, g = 0, d = 0){
  eff_depreciation <- k*(n + g + d)
  return(eff_depreciation)
}

country_tbl <- country_tbl %>% 
  mutate(savings_curve = find.saving.point(k = k,
                                           alpha = alpha,
                                           s = s),
         depreciation_curve = find.effective.depreciation(k = k,
                                                          n = n_g_d),
         k_star = (s / n_g_d) ^ (1 / (1 - alpha)),
         s_star = find.saving.point(k = k_star,
                                    alpha = alpha,
                                    s = s))
country_tbl %>% 
  sample_n(5) %>% 
  kable(digits = 2)
kcountry1985sn_g_dclubmodel_typealphabetasavings_curvedepreciation_curvek_stars_star
16.5Hong Kong133720.200.08intermediate_no_oilAugmented0.300.250.461.323.680.29
8.0Sri Lanka24820.150.07intermediate_no_oilAugmented0.300.250.280.592.690.20
13.0Hong Kong133720.200.08intermediate_no_oilAugmented0.300.250.431.043.680.29
12.5MaltaNA0.31NAoilAugmented0.99-1.803.77NANANA
17.5Liberia9440.220.08developing_no_oilClassic0.36NA0.601.404.690.37

We now have all the information we need to plot each countries’ curves. To demonstrate the difference between subsets we use facet_wrap(~club). The Crosstalk library lets us highlight individual countries within a plot and compare across multiple plots throughout the post.

First, the classic Solow model savings curves:

[Edit: Plotly is having some issues with GitHub Pages’ Liquid tags currently so I’m using an iframe workaround that doesn’t allow for crosstalk’s filter functionality so the plots below are quite messy.]

library(crosstalk)
plot_dat <- country_tbl %>% 
         filter(model_type == "Classic") %>% 
         filter(club != "oil")

plot_dat <- SharedData$new(plot_dat, ~country, group = "Choose a Country:")

p <- ggplot(plot_dat,
            aes(k,
                savings_curve,
                color = country,
                text = country)) +
  geom_line() +
  theme_minimal() +
  facet_wrap(~club)
q <- p %>% 
  ggplotly(tooltip = "text") %>%   
  highlight(on = "plotly_click",
            persistent = TRUE,
            selectize = TRUE,
            defaultValues = c("United Kingdom",
                              "Angola",
                              "Brazil"),
            opacity = 0.1) %>% 
  layout(showlegend = FALSE)
q

The differences between the rich, developing and intermediate savings curves demonstrate visually that the Solow model predicts conditional rather than unconditional convergence. The fact that intermediate countries’ savings curves appear to be higher than the rich countries, whilst puzzling, is coherent with the Solow model.

Now, we add depreciation curves as well as markers for \(k^*\) and \(s^*\):

## Now adding depreciation curves and markers
p_dep <- p +
  geom_line(aes(k, depreciation_curve)) +
  geom_point(aes(k_star, s_star))

q <- p_dep %>% 
  ggplotly(tooltip = "text") %>%   
  highlight(on = "plotly_click",
            persistent = TRUE,
            selectize = TRUE,
            defaultValues = c("United Kingdom",
                              "Angola",
                              "Brazil"),
            opacity = 0.1) %>% 
  layout(showlegend = FALSE)
q

The above plot is a little messy, so here we just show the curves’ intersection:

plot_dat <- country_tbl %>% 
         filter(model_type == "Classic") %>% 
         filter(club != "oil") %>% 
         select(country, k_star, s_star, club) %>% 
         distinct()

plot_dat <- SharedData$new(plot_dat, ~country, group = "Choose a Country:")
p <- ggplot(plot_dat,
            aes(k_star, s_star, color = country, text = country)) +
  geom_point() +
  facet_wrap(~club) +
  theme_minimal()
q <- p %>% 
  ggplotly(tooltip = "text") %>%   
  highlight(on = "plotly_click",
            persistent = TRUE,
            selectize = TRUE,
            defaultValues = c("United Kingdom",
                              "Angola",
                              "Brazil"),
            opacity = 0.1) %>% 
  layout(showlegend = FALSE)
q

Plotting augmented saving curves

Now we can do the same but with the Augmented Solow model:

Comparing the two together:

This plot supports MRW’s conclusion that human capital can help explain differences in income. The addition of school in the regression equation and therefore \(\beta\) in the Augmented Solow model now means that rich countries’ savings curves appear above the intermediate countries as we’d most likely expect. Furthermore, it seems like the poor countries have gotten poorer.

Conclusion

We’ve plotted the savings and depreciation curves of the Solow and Augmented Solow model. Whilst the plots are interactive to a degree, thanks to Plotly and Crosstalk, it’d be interesting to allow users to input their own parameters for \(\alpha\), \(\beta\), \(s\), \((n + g + d)\) - unfortunately this would require Shiny which I’m reluctant to use given GitHub Pages’ static hosting.

Again our visualisations and tweaks of MRW’s original paper broadly reflect their own findings. Part I of this in depth MRW marathon can be found here, followed by Part II and Part III.