From 1e0ba7b4ff9b71f3e8aaaee98acbd115ec66592b Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Thu, 24 Oct 2024 16:42:48 -0700 Subject: [PATCH] update make_dsem_ram to use additive-kronecker --- NAMESPACE | 3 ++ R/dsem.R | 2 +- R/make_dsem_ram.R | 52 +++++++++++++++++++++++++--------- tests/testthat/test-platform.R | 16 +++++------ 4 files changed, 50 insertions(+), 23 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index ca85fe8..eda8d97 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,7 +22,10 @@ export(make_dsem_ram) export(parse_path) export(stepwise_selection) importFrom(Matrix,Cholesky) +importFrom(Matrix,drop0) +importFrom(Matrix,mat2triplet) importFrom(Matrix,solve) +importFrom(Matrix,sparseMatrix) importFrom(TMB,MakeADFun) importFrom(TMB,compile) importFrom(TMB,dynlib) diff --git a/R/dsem.R b/R/dsem.R index a90aa8e..e6dd227 100644 --- a/R/dsem.R +++ b/R/dsem.R @@ -25,7 +25,7 @@ #' #' @importFrom TMB compile dynlib MakeADFun sdreport summary.sdreport #' @importFrom stats AIC sd .preformat.ts na.omit nlminb optimHess pnorm rnorm simulate time tsp<- -#' @importFrom Matrix solve Cholesky +#' @importFrom Matrix solve Cholesky sparseMatrix mat2triplet drop0 #' @importFrom sem sem #' @importFrom igraph plot.igraph graph_from_data_frame with_sugiyama layout_ #' @importFrom ggraph ggraph geom_edge_arc create_layout rectangle geom_node_text theme_graph diff --git a/R/make_dsem_ram.R b/R/make_dsem_ram.R index dbbe6e4..e31250e 100644 --- a/R/make_dsem_ram.R +++ b/R/make_dsem_ram.R @@ -332,22 +332,46 @@ function( sem, } # Loop through paths + G_kk = P_kk = sparseMatrix( i=1, j=1, x=0, dims=rep(length(variables)*length(times),2) ) # Make with a zero + #P_kk = new("dgCMatrix") + #P_kk = Matrix() + #P_kk@Dim <- as.integer(rep(length(variables)*length(times),2)) + #P_kk@p = integer(length(variables)*length(times)+1L) + #G_kk = P_kk for( i in seq_len(nrow(model)) ){ - for( t in seq_along(times) ){ lag = as.numeric(model[i,2]) - par.no = par.nos[i] - # Get index for "from" - from = c( times[t], model[i,'first'] ) - from_index = match_row( Q_names, from ) - from_index = ifelse( length(from_index)==0, NA, from_index ) - # Get index for "to" - to = c( times[t+lag], model[i,'second'] ) - to_index = match_row( Q_names, to ) - to_index = ifelse( length(to_index)==0, NA, to_index ) - ram_new = data.frame( "heads"=abs(as.numeric(model[i,'direction'])), "to"=to_index, "from"=from_index, "parameter"=par.no, "start"=startvalues[i] ) - ram = rbind( ram, ram_new ) - }} - rownames(ram) = NULL + L_tt = sparseMatrix( i = seq(lag+1,length(times)), + j = seq(1,length(times)-lag), + x = 1, + dims = rep(length(times),2) ) + P_jj = sparseMatrix( i = match(model[i,'second'],variables), + j = match(model[i,'first'],variables), + x = 1, + dims = rep(length(variables),2) ) + tmp_kk = kronecker(P_jj, L_tt) + if(abs(as.numeric(model[i,'direction']))==1){ + P_kk = P_kk + tmp_kk * par.nos[i] + }else{ + G_kk = G_kk + tmp_kk * par.nos[i] + } + #for( t in seq_along(times) ){ + # # Get index for "from" + # from = c( times[t], model[i,'first'] ) + # from_index = match_row( Q_names, from ) + # from_index = ifelse( length(from_index)==0, NA, from_index ) + # # Get index for "to" + # to = c( times[t+lag], model[i,'second'] ) + # to_index = match_row( Q_names, to ) + # to_index = ifelse( length(to_index)==0, NA, to_index ) + # ram_new = data.frame( "heads"=abs(as.numeric(model[i,'direction'])), "to"=to_index, "from"=from_index, "parameter"=par.no, "start"=startvalues[i] ) + # ram = rbind( ram, ram_new ) + #} + } + #rownames(ram) = NULL + ram = rbind( cbind(1, sapply(mat2triplet(drop0(P_kk)),cbind)), + cbind(2, sapply(mat2triplet(drop0(G_kk)),cbind)) ) + ram = data.frame( ram, startvalues[ram[,4]] ) + colnames(ram) = c("heads", "to", "from", "parameter", "start") # if( isTRUE(remove_na) ){ diff --git a/tests/testthat/test-platform.R b/tests/testthat/test-platform.R index 3987d6c..5460773 100644 --- a/tests/testthat/test-platform.R +++ b/tests/testthat/test-platform.R @@ -3,23 +3,23 @@ test_that("dsem example is working ", { #skip_on_ci() sem = " Profits -> Consumption, 0, a2 - Profits -> Consumption, -1, a3 + Profits -> Consumption, 1, a3 Priv_wage -> Consumption, 0, a4 Gov_wage -> Consumption, 0, a4 Consumption <-> Consumption, 0, v1 - Consumption -> Consumption, -1, ar1 - Consumption -> Consumption, -2, ar2 + Consumption -> Consumption, 1, ar1 + Consumption -> Consumption, 2, ar2 Profits -> Investment, 0, b2 - Profits -> Investment, -1, b3 - Capital_stock -> Investment, -1, b4 + Profits -> Investment, 1, b3 + Capital_stock -> Investment, 1, b4 Investment <-> Investment, 0, v2 neg_Gov_wage <-> neg_Gov_wage, 0, v3 GNP -> Priv_wage, 0, c2 Taxes -> Priv_wage, 0, c2 neg_Gov_wage -> Priv_wage, 0, c2 - GNP -> Priv_wage, -1, c3 - Taxes -> Priv_wage, -1, c3 - neg_Gov_wage -> Priv_wage, -1, c3 + GNP -> Priv_wage, 1, c3 + Taxes -> Priv_wage, 1, c3 + neg_Gov_wage -> Priv_wage, 1, c3 Time -> Priv_wage, 0, c4 Priv_wage <-> Priv_wage, 0, v4 GNP <-> GNP, 0, v5