Approximation using multi-steps solutions

[paulofelipe]

First of all, thanks for making this package available. It is a great initiative!

Unfortunately, I was unable to see the presentation at the GTAP conference. However, yesterday I decided to test the package. And it's impressive!

As you show in the paper, there are differences in the results of the solution in steps and that of the gempack (gragg 2-4-6). In my tests, there is also a difference for the case of a single solution that uses steps and the gempack solution (euler 1 solution). I believe there is something going on internally in the solution function, but I still haven't been able to identify exactly.

To test, I created two functions that use the package's one step solution. Within these functions, sub-shocks are created and the accumulated results are generated. In this way, I was able to approximate the results of the GTAP. I believe it would be interesting to test these functions in your example, but I do not have access to the GTAP 10 base. I think that later it will be possible to implement a function that exactly replicates the Gragg method. For now, the function is closer to the Euler method.

Anyway, follow my example (which needs to be improved) below. It uses the example "BOOK3x3" available in runGTAP.

Thanks!

library(HARr)
library(tabloToR)
library(tidyverse)

# Functions to solutions in steps  ----------------------------------------------

solve_model_steps <- function(model, shocks, steps){

  model <- model
  data <- model$data
   
  sub_shocks <- ifelse(
    names(shocks) %in% model$changeVariables,
    shocks/iter,
    ((1 + shocks/100)^(1/steps) -  1) * 100
  )

  names(sub_shocks) <- names(shocks)

  model$setShocks(sub_shocks)

  all_variables <- unique(str_extract(model$data$variables, "^.*?(?=\\[)"))
  basic_change_variables <- model$basicChangeVariables
  sol <- list()

  for (j in 1:steps) {
    model$solveModel(iter = 1)
  
    if (j == 1) {
      for (i in all_variables) {
        sol[[i]] <- model$data[[i]]
        model$data[[i]] <- model$data[[i]] - model$data[[i]]
      }
    }
  
    if (j > 1) {
      for (i in all_variables) {
        if (i %in% basic_change_variables) {
          sol[[i]] <- sol[[i]] + model$data[[i]]
        } else {
          sol[[i]] <- ((1 + sol[[i]] / 100) * (1 + model$data[[i]] / 100) - 1) * 100
        }
        model$data[[i]] <- model$data[[i]] - model$data[[i]]
      }
    }
  }

  model$data <- data
  
  for (i in all_variables) {
    model$data[[i]] <- sol[[i]]
  }

  model$data <- model$generateUpdates(model$data)

  return(model)

}

# Steps should be 1N,  2N, 4N
solve_model_multisteps <- function(model, data, shocks, steps = c(1, 2, 4)){

  all_variables <- unique(str_extract(model$data$variables, "^.*?(?=\\[)"))
  model$loadData(data)
  sol1 <- solve_model_steps(model, shocks, steps[1])
  data1 <- sol1$data
  model$loadData(data)
  sol2 <- solve_model_steps(model, shocks, steps[2])
  data2 <- sol2$data
  model$loadData(data)
  sol3 <- solve_model_steps(model, shocks, steps[3])
  data3 <- sol1$data

  #model <- sol1
  model$loadData(data)
  for (i in all_variables) {
    model$data[[i]] <- (8 * data3[[i]] - 6 * data2[[i]] + data1[[i]])/3
  }

  model$data <- model$generateUpdates(model$data)

  return(model)

}

# Create the model -----------------------------------------------------------

GTAP <- GEModel$new()
GTAP$loadTablo("tablo/gtap.tab")
#> Carregando pacotes exigidos: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack

# Load data ------------------------------------------------------------------
data = list(
  GTAPSETS = read_har("data/sets.har"),
  GTAPPARM = read_har("data/default.prm"),
  GTAPDATA = read_har("data/basedata.har")
)
#> Scanning records for headers

#> Found 8 headers
#> Sorting out records within headers
#> Processing first and second records within headers
#> Processing character headers
#> Processing integer headers
#> Processing real headers
#> Scanning records for headers

#> Found 14 headers
#> Sorting out records within headers
#> Processing first and second records within headers
#> Processing character headers
#> Processing integer headers
#> Processing real headers
#> Scanning records for headers

#> Found 29 headers
#> Sorting out records within headers
#> Processing first and second records within headers
#> Processing character headers
#> Processing integer headers
#> Processing real headers

GTAP$loadData(data)

# Closure --------------------------------------------------------------------

allVariables <- GTAP$data$variables

exogenousVariables <- c(
  "afall", "afcom", "afeall", "afecom", "afereg", "afesec", "afreg", "afsec",
  "ams", "aoall", "aoreg", "aosec", "atall", "atd", "atf", "atm", "ats", "au",
  "avaall", "avareg", "avasec", "cgdslack", "dpgov", "dppriv", "dpsave",
  "endwslack", "incomeslack", "pfactwld", "pop", "profitslack", "psaveslack",
  "tf", "tfd", "tfm", "tgd", "tgm", "tm", "tms", "to", "tpd", "tpm", "tp",
  "tradslack", "tx", "txs"
)

exogenousModelVariables <- allVariables[(Reduce(function(a, f) {
    c(a, grep(sprintf("^%s\\[", f), allVariables))
  }, exogenousVariables, c())
)]

for (r in GTAP$data$REG) {
  for (e in GTAP$data$ENDW_COMM) {
    exogenousModelVariables <- c(exogenousModelVariables, sprintf('qo["%s","%s"]', e, r))
  }
}

# Shocks ----------------------------------------------------------------------

shocks <- array(0, dim = length(exogenousModelVariables), dimnames = list(exogenousModelVariables))
shocks['tms["food","USA","EU"]'] <- -10

# Solve using euler multisteps -------------------------------------------------

GTAP <- solve_model_multisteps(GTAP, data, shocks, steps = c(4, 8, 16))


round(GTAP$data$qgdp, 6)
#>       USA        EU       ROW 
#> -0.000164  0.011983 -0.001232
round(GTAP$data$EV, 2)
#>     USA      EU     ROW 
#>  886.55  216.48 -439.58

# GTAP Solution Gragg 2-4-6
# EV
# 0175  EV
# 1 USA 886,60
# 2 EU  216,24
# 3 ROW -439,59
# Total 663,25

# qgdp
# 0136  qgdp
# 1 USA -0,000163
# 2 EU  0,011981
# 3 ROW -0,001229

^{Created on 2020-06-20 by the reprex package (v0.3.0)}