Full incidence matrix about bimets HOT 3 CLOSED

rm(list=ls(all.names = T))

library(bimets)
library(igraph)
library(extrafont)
library(hrbrthemes)
library(ggalt)
library(ggplot2)
library(plotly)

#define model
klein1_def <- "
MODEL

COMMENT> Consumption
BEHAVIORAL> cn
TSRANGE 1921 1 1941 1
EQ> cn = a1 + a2*p + a3*TSLAG(p,1) + a4*(w1+w2)
COEFF> a1 a2 a3 a4

COMMENT> Investment
BEHAVIORAL> i
TSRANGE 1921 1 1941 1
EQ> i = b1 + b2*p + b3*TSLAG(p,1) + b4*TSLAG(k,1)
COEFF> b1 b2 b3 b4

COMMENT> Demand for Labor
BEHAVIORAL> w1
TSRANGE 1921 1 1941 1
EQ> w1 = c1 + c2*(y+t-w2) + c3*TSLAG(y+t-w2,1)+c4*time
COEFF> c1 c2 c3 c4

COMMENT> Gross National Product
IDENTITY> y
EQ> y = cn + i + g - t 

COMMENT> Profits
IDENTITY> p
EQ> p = y - (w1+w2)

COMMENT> Capital Stock 
IDENTITY> k
EQ> k = TSLAG(k,1) + i
END
"

#load model
myModel <- LOAD_MODEL(modelText = klein1_def)

#define vertices, i.e. vendog + vexog
vertices_names  <- c(myModel$vendog,myModel$vexog)
vertices_type   <- c(rep("en",length(myModel$vendog)),rep("ex",length(myModel$vexog)))

vertices  <- data.frame(name=vertices_names,
                        type=vertices_type)

#define relations between vertices
relations_from  <- c()
relations_to    <- c()

for (idx in names(myModel$behaviorals))
{
  #get eq components names
  localIncidenceVars <- myModel$behaviorals[[idx]]$eqComponentsNames
  #remove dependant 
  localIncidenceVars <- setdiff(localIncidenceVars,idx)
  
  for (from in localIncidenceVars)
  {
    relations_from  <- c(relations_from,from)
    relations_to    <- c(relations_to,idx)
  }
}

for (idx in names(myModel$identities))
{
  localIncidenceVars <- myModel$identities[[idx]]$eqComponentsNames
  localIncidenceVars <- setdiff(localIncidenceVars,idx)
  
  for (from in localIncidenceVars)
  {
    relations_from  <- c(relations_from,from)
    relations_to    <- c(relations_to,idx)
  }
}
 
relations <- data.frame(from=relations_from,
                        to=relations_to)


#add attribute to graph, TRUE if from-node is an endogenous
relations$fromIsEndogenous <- FALSE
relations$fromIsEndogenous[which(relations$from %in% vertices$name[which(vertices$type=='en')])] <- TRUE

#create graph
g <- graph.data.frame(relations, directed=TRUE, vertices=vertices)

#define edge colors, dark blue if circular
edge_colors <- rep("pink",length(relations$from))
edge_colors[which(E(g)$fromIsEndogenous==TRUE)] <- "skyblue"

#plot stuff
set.seed(18)
plot(g,
     main="Klein model full incidence graph",
     edge.arrow.size=0.3, 
     edge.arrow.width=3,
     edge.width=c(1.3, 2.1)[1+(E(g)$fromIsEndogenous==TRUE)],
     edge.color=edge_colors,
     vertex.color=c("pink", "skyblue")[1+(V(g)$type=="en")],
     vertex.size=c(35,30)[1+(V(g)$type=="en")], 
     vertex.frame.color="#ffffff",
     vertex.label.color=c("#662222", "#222266")[1+(V(g)$type=="en")],
     vertex.label.cex=0.8, 
     vertex.label.family="Arial Black", 
     edge.curved=-0.2,
     layout=layout.davidson.harel
)