Added constant Elasticity code

CobbDouglas/LICENSE

@@ -7,7 +7,7 @@ of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights
 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
+furnished to do so, subject to the following conditions-
 
 The above copyright notice and this permission notice shall be included in all
 copies or substantial portions of the Software.

CobbDouglas/README.md

@@ -1,22 +1,22 @@
 ## README
 
 ## About 
-This is some R code for the paper [Labor, Capital, and the Optimal Growth ofSocial Movements](https://nunosempere.github.io/ea/MovementBuildingForUtilityMaximizers.pdf)
+This is some R code for the paper [Labor, Capital, and the Optimal Growth ofSocial Movements](https-//nunosempere.github.io/ea/MovementBuildingForUtilityMaximizers.pdf)
 
 ## Structure of the code.
 
-1. `Variables.R`: Contains variables
-2. `Transition dynamics.R`: If we know our system at time t, generate an approximation of our system at time t+stepsize. 
-3. `Forward Shooting.R`: Carries out the forward shooting. 
-4. `ReverseShooting.R`: Nonfunctional. Ignore. 
-5. `Plotting.R`: Generates graphs of the results
+1. `Variables.R`- Contains variables
+2. `Transition dynamics.R`- If we know our system at time t, generate an approximation of our system at time t+stepsize. 
+3. `Forward Shooting.R`- Carries out the forward shooting. 
+4. `ReverseShooting.R`- Nonfunctional. Ignore. 
+5. `Plotting.R`- Generates graphs of the results
 
 ## How to run
 
 Open files 1.,2.,3., and 5. in an IDE for R, like RStudio. Run them in order. For `Plotting.R`, add or create the directory in which you want the graphs to be generated, making sure it has suitable permissions. 
 
 ## Gotchas
-Using a very small stepsize runs into floating point errors. Consider a stylized example:
+Using a very small stepsize runs into floating point errors. Consider a stylized example-
 
 ```r
 options(digits=22)
@@ -28,7 +28,7 @@ stepsize <- 10^(-3)
 ## Example 1
 x <- pi*1e+60
 print(x)
-for(i in c(1:numsteps)){
+for(i in c(1-numsteps)){
   x <- x+dx*stepsize
 }
 x == pi*1e+60
@@ -45,20 +45,20 @@ The two examples should give the same results, but don't.
 
 ### Why reverse shooting doesn't work
 
-Consider this example:
+Consider this example-
 
 ```r
 ## Stylized forward shooting
 
 x <- 0
-for(i in c(1:200)){
+for(i in c(1-200)){
   x <- x + 7^i
 }
 
 ## Stylized reverse shooting
 
 y <- x
-for(i in c(200:1)){
+for(i in c(200-1)){
   y <- y - 7^i
 }
 print(y)

CobbDouglas/code/ForwardShooting.R

@@ -59,5 +59,5 @@ a3_growth = (a3_array_forward_shooting[l]-a3_array_forward_shooting[l-1])/a3_arr
 a3_growth
 
 a1_array_forward_shooting[l]/x1_array_forward_shooting[l]
-plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100):l])
+plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100)-l])
 

CobbDouglas/code/Plotting.R

@@ -42,7 +42,7 @@ saveplot = function(imagename){
   if(saveplots){
     ggsave(paste(imagenumbercounter,"_", imagename, "_",shootingtype, "shooting",".png", sep =""), units="in", width=width, height=height)  
     imagenumbercounter <<- imagenumbercounter+1
-    ## https://stackoverflow.com/questions/1236620/global-variables-in-r
+    ## https-//stackoverflow.com/questions/1236620/global-variables-in-r
   }
 }
 options(digits=1) ## Just for display
@@ -77,7 +77,7 @@ title_text = "Evolution of state variables"
 )
 saveplot("StateVariablesX1X3")
 
-## x1 and x3: log plot
+## x1 and x3- log plot
 title_text = "Evolution of state variables"
 (ggplot(data = xs, aes(x = times, y= values, color = var)) +
     geom_line(size = 0.5)
@@ -408,7 +408,7 @@ title_text = "Evolution of spending"
 )
 saveplot("SpendingA1A3")
 
-## alpha1 and alpha3: log plot
+## alpha1 and alpha3- log plot
 
 alphas <- list()
 alphas$values <- c(a1_array_plotting, a3_array_plotting)

CobbDouglas/code/ReverseShooting.R

@@ -35,7 +35,7 @@ for(t in times_reverse_shooting){
   x1_t = (x1_t - dx1(t-stepsize)*stepsize)/(1+r1_stepsize)
   
   ## But then I can calulate s1, s3 again, and do another iteration.
-  # for(i in c(1:10)){
+  # for(i in c(1-10)){
   #   # using the new x3_t
   #   s3_t = s3(t_previous)
   #   x3_t = dx3_reverse_shooting(t, guessmin, guessmax, x3_target) 

CobbDouglas/code/Variables.R

@@ -5,7 +5,7 @@ r_1 = 0.06
 r_3 = -0.02
 γ_1 = 0.03
 γ_3 = 0.01
-w_3 = 5000 ## Mean on the EA Survey 2018: ~7000
+w_3 = 5000 ## Mean on the EA Survey 2018- ~7000
 β_3 = 1 ## 1 ## 0.5 Corresponds to roughly 5 people convincing 10 other people per year on a 50k year budget
   ## Dramatic if beta = 0.5, w_3 = 1000
 λ_1 = 0.5

CobbDouglas/plots/10000years/.Rhistory

@@ -61,7 +61,7 @@ a1_growth
 a3_growth = (a3_array_forward_shooting[l]-a3_array_forward_shooting[l-1])/a3_array_forward_shooting[l-1]/stepsize
 a3_growth
 a1_array_forward_shooting[l]/x1_array_forward_shooting[l]
-plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100):l])
+plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100)-l])
 # Plotting
 ## install.packages("tidyverse") <- Not totally necessary.
 ## install.packages("ggplot2")
@@ -100,7 +100,7 @@ saveplot = function(imagename){
 if(saveplots){
 ggsave(paste(imagenumbercounter,"_", imagename, "_",shootingtype, "shooting",".png", sep =""), units="in", width=width, height=height)
 imagenumbercounter <<- imagenumbercounter+1
-## https://stackoverflow.com/questions/1236620/global-variables-in-r
+## https-//stackoverflow.com/questions/1236620/global-variables-in-r
 }
 }
 options(digits=1) ## Just for display
@@ -135,7 +135,7 @@ labels=c("Capital", "Labor")
 + scale_y_continuous(breaks = seq(min(x1_array_plotting), max(x1_array_plotting), length.out=5))
 )
 saveplot("StateVariablesX1X3")
-## x1 and x3: log plot
+## x1 and x3- log plot
 title_text = "Evolution of state variables"
 (ggplot(data = xs, aes(x = times, y= values, color = var)) +
 geom_line(size = 0.5)
@@ -478,7 +478,7 @@ labels=c("Direct spending","Spending on movement building")
 )
 )
 saveplot("SpendingA1A3")
-## alpha1 and alpha3: log plot
+## alpha1 and alpha3- log plot
 alphas <- list()
 alphas$values <- c(a1_array_plotting, a3_array_plotting)
 alphas$fraction <- c(rep("a1", length(s1_array_plotting)), rep("a2", length(s3_array_plotting)))

CobbDouglas/plots/1000years/.Rhistory

@@ -8,7 +8,7 @@ a1_growth
 a3_growth = (a3_array_forward_shooting[l]-a3_array_forward_shooting[l-1])/a3_array_forward_shooting[l-1]/stepsize
 a3_growth
 a1_array_forward_shooting[l]/x1_array_forward_shooting[l]
-plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100):l])
+plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100)-l])
 # Forward shooting
 options(digits=7)
 ## Evolution
@@ -48,7 +48,7 @@ r_1 = 0.06
 r_3 = -0.05
 γ_1 = 0.03
 γ_3 = 0.01
-w_3 = 5000 ## Mean on the EA Survey 2018: ~7000
+w_3 = 5000 ## Mean on the EA Survey 2018- ~7000
 β_3 = 1 ## 1 ## 0.5 Corresponds to roughly 5 people convincing 10 other people per year on a 50k year budget
 ## Dramatic if beta = 0.5, w_3 = 1000
 λ_1 = 0.5
@@ -120,7 +120,7 @@ a1_growth
 a3_growth = (a3_array_forward_shooting[l]-a3_array_forward_shooting[l-1])/a3_array_forward_shooting[l-1]/stepsize
 a3_growth
 a1_array_forward_shooting[l]/x1_array_forward_shooting[l]
-plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100):l])
+plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100)-l])
 # Plotting
 ## install.packages("tidyverse") <- Not totally necessary.
 ## install.packages("ggplot2")
@@ -158,7 +158,7 @@ saveplot = function(imagename){
 if(saveplots){
 ggsave(paste(imagenumbercounter,"_", imagename, "_",shootingtype, "shooting",".png", sep =""), units="in", width=width, height=height)
 imagenumbercounter <<- imagenumbercounter+1
-## https://stackoverflow.com/questions/1236620/global-variables-in-r
+## https-//stackoverflow.com/questions/1236620/global-variables-in-r
 }
 }
 options(digits=1) ## Just for display
@@ -191,7 +191,7 @@ labels=c("Capital", "Labor")
 )
 + scale_y_continuous(breaks = seq(min(x1_array_plotting), max(x1_array_plotting), length.out=5))
 )
-## x1 and x3: log plot
+## x1 and x3- log plot
 title_text = "Evolution of state variables"
 (ggplot(data = xs, aes(x = times, y= values, color = var)) +
 geom_line(size = 0.5)
@@ -282,7 +282,7 @@ r_1 = 0.06
 r_3 = -0.02
 γ_1 = 0.03
 γ_3 = 0.01
-w_3 = 5000 ## Mean on the EA Survey 2018: ~7000
+w_3 = 5000 ## Mean on the EA Survey 2018- ~7000
 β_3 = 1 ## 1 ## 0.5 Corresponds to roughly 5 people convincing 10 other people per year on a 50k year budget
 ## Dramatic if beta = 0.5, w_3 = 1000
 λ_1 = 0.5
@@ -354,7 +354,7 @@ a1_growth
 a3_growth = (a3_array_forward_shooting[l]-a3_array_forward_shooting[l-1])/a3_array_forward_shooting[l-1]/stepsize
 a3_growth
 a1_array_forward_shooting[l]/x1_array_forward_shooting[l]
-plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100):l])
+plot((a1_array_forward_shooting/x1_array_forward_shooting)[(l-100)-l])
 # Plotting
 ## install.packages("tidyverse") <- Not totally necessary.
 ## install.packages("ggplot2")
@@ -392,7 +392,7 @@ saveplot = function(imagename){
 if(saveplots){
 ggsave(paste(imagenumbercounter,"_", imagename, "_",shootingtype, "shooting",".png", sep =""), units="in", width=width, height=height)
 imagenumbercounter <<- imagenumbercounter+1
-## https://stackoverflow.com/questions/1236620/global-variables-in-r
+## https-//stackoverflow.com/questions/1236620/global-variables-in-r
 }
 }
 options(digits=1) ## Just for display
@@ -426,7 +426,7 @@ labels=c("Capital", "Labor")
 + scale_y_continuous(breaks = seq(min(x1_array_plotting), max(x1_array_plotting), length.out=5))
 )
 saveplot("StateVariablesX1X3")
-## x1 and x3: log plot
+## x1 and x3- log plot
 title_text = "Evolution of state variables"
 (ggplot(data = xs, aes(x = times, y= values, color = var)) +
 geom_line(size = 0.5)

ConstantElasticity/1-Variables.R

@@ -0,0 +1,44 @@
+## Variables
+η = 1.1
+ρ = -0.5
+δ = 0.005
+r_1 = 0.06
+r_2 = -0.02
+### r_1_paper = 0.06 ## yearly
+### r_1 = 1+log(1+r_1_paper) ## instantaneous
+### r_2_paper = -0.02 ## yearly
+### r_2 = 1+log(1+r_2_paper) ## instantaneous
+γ_1 = 0.03
+γ_2 = 0.01
+w_2 = 5000 
+β_2 = 1 
+λ_1 = 0.5
+λ_2 = 0.5
+δ_2 = 0.44
+q=0.5
+
+## Initial conditions. Correspond to "year 0". 
+K_init = 10^13 ## 10^13 to afford movement building. 10^14 to have a reasonable amount of direct investment as well. This is close to the net discounted value of US GDP.
+L_init = 10^4 
+
+## Integration constant
+c1_forward_shooting = 10^(-8) # As low as possible
+
+## Stepsize
+stepsize = 0.1
+r1_stepsize = ((1+r_1)^stepsize)-1 
+r2_stepsize = ((1+r_2)^stepsize)-1 
+
+## Notes time it takes to run the simulations
+### stepsize = 0.1 => seconds (7s). 
+### stepsize = 0.01 => minutes (3 mins).
+
+## Knife-edge constant
+δ  = r_1 + η*( (γ_1/(ρ-1)) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2)) )
+knife_edge_constant = (γ_1/(ρ-1)) + ((r_1 - δ)/η) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2))
+
+## Interval
+first = 0
+last = 1000
+times_forward_shooting = seq(from=first, to=last, by=stepsize)
+times_reverse_shooting = seq(from=last, to=first, by=-stepsize)

ConstantElasticity/2-Scenarios.R

@@ -0,0 +1,167 @@
+# ========================================================== #
+#                      scenario1                             #
+# ========================================================== #
+# Description: 
+
+## Variables
+η = 1.1
+ρ = -0.5
+δ = 0.005
+r_1 = 0.06
+r_2 = -0.02
+### r_1_paper = 0.06 ## yearly
+### r_1 = 1+log(1+r_1_paper) ## instantaneous
+### r_2_paper = -0.02 ## yearly
+### r_2 = 1+log(1+r_2_paper) ## instantaneous
+γ_1 = 0.03
+γ_2 = 0.01
+w_2 = 5000 
+β_2 = 1 
+λ_1 = 0.5
+λ_2 = 0.5
+δ_2 = 0.44
+q=0.5
+
+## Initial conditions. Correspond to "year 0". 
+K_init = 10^13
+L_init = 10^4 
+
+## Integration constant
+c1_forward_shooting = 10^(-8) # 10^10 ## 2*10^(-6)
+
+## Stepsize
+stepsize = 0.1
+r1_stepsize = ((1+r_1)^stepsize)-1 
+r2_stepsize = ((1+r_2)^stepsize)-1 
+
+## Notes time it takes to run the simulations
+### stepsize = 0.1 => seconds (7s). 
+### stepsize = 0.01 => minutes (3 mins).
+
+## Knife-edge constant
+knife_edge_constant = (γ_1/(ρ-1)) + ((r_1 - δ)/η) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2))
+
+## Interval
+first = 0
+last = 1000
+times_forward_shooting = seq(from=first, to=last, by=stepsize)
+times_reverse_shooting = seq(from=last, to=first, by=-stepsize)
+
+# ========================================================== #
+#                      scenario2                             #
+# ========================================================== #
+# Description: Knife edge constant = 0
+
+## Variables
+η = 1.1
+ρ = -0.5
+### δ = 0.005
+r_1 = 0.06
+r_2 = -0.02
+### r_1_paper = 0.06 ## yearly
+### r_1 = 1+log(1+r_1_paper) ## instantaneous
+### r_2_paper = -0.02 ## yearly
+### r_2 = 1+log(1+r_2_paper) ## instantaneous
+γ_1 = 0.03
+γ_2 = 0.01
+w_2 = 5000 
+β_2 = 1 
+λ_1 = 0.5
+λ_2 = 0.5
+δ_2 = 0.44
+q=0.5
+δ  = r_1 + η*( (γ_1/(ρ-1)) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2)) )
+
+## Initial conditions. Correspond to "year 0". 
+K_init = 10^13 ## 10^13 to afford movement building. 10^14 to have a reasonable amount of direct investment as well. This is close to the net discounted value of US GDP.
+L_init = 10^4 
+
+## Integration constant
+c1_forward_shooting = 10^(-8) # 10^10 ## 2*10^(-6)
+
+## Stepsize
+stepsize = 0.1
+r1_stepsize = ((1+r_1)^stepsize)-1 
+r2_stepsize = ((1+r_2)^stepsize)-1 
+
+## Notes time it takes to run the simulations
+### stepsize = 0.1 => seconds (7s). 
+### stepsize = 0.01 => minutes (3 mins).
+
+## Knife-edge constant
+knife_edge_constant = (γ_1/(ρ-1)) + ((r_1 - δ)/η) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2))
+
+## Interval
+first = 0
+last = 1000
+times_forward_shooting = seq(from=first, to=last, by=stepsize)
+times_reverse_shooting = seq(from=last, to=first, by=-stepsize)
+
+
+# ========================================================== #
+#                      scenario3                             #
+# ========================================================== #
+# Description: Knife edge constant < 0
+
+## Variables
+η = 1.1
+ρ = -0.5
+δ = 0.00578572
+r_1 = 0.06
+r_2 = -0.02
+### r_1_paper = 0.06 ## yearly
+### r_1 = 1+log(1+r_1_paper) ## instantaneous
+### r_2_paper = -0.02 ## yearly
+### r_2 = 1+log(1+r_2_paper) ## instantaneous
+γ_1 = 0.03
+γ_2 = 0.01
+w_2 = 5000 
+β_2 = 1 
+λ_1 = 0.5
+λ_2 = 0.5
+δ_2 = 0.44
+q=0.5
+
+## Initial conditions. Correspond to "year 0". 
+K_init = 10^13 ## 10^13 to afford movement building. 10^14 to have a reasonable amount of direct investment as well. This is close to the net discounted value of US GDP.
+L_init = 10^4 
+
+## Integration constant
+c1_forward_shooting = 10^(-8) # 10^10 ## 2*10^(-6)
+
+## Stepsize
+stepsize = 0.1
+r1_stepsize = ((1+r_1)^stepsize)-1 
+r2_stepsize = ((1+r_2)^stepsize)-1 
+
+## Notes time it takes to run the simulations
+### stepsize = 0.1 => seconds (7s). 
+### stepsize = 0.01 => minutes (3 mins).
+
+## Knife-edge constant
+knife_edge_constant = (γ_1/(ρ-1)) + ((r_1 - δ)/η) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2))
+knife_edge_constant
+## Interval
+first = 0
+last = 1000
+times_forward_shooting = seq(from=first, to=last, by=stepsize)
+times_reverse_shooting = seq(from=last, to=first, by=-stepsize)
+
+# ========================================================== #
+#                     scenarios 4 to 14                      #
+# ========================================================== #
+
+δ_array = seq(from=0.0044, to=0.0064, by=0.0002)
+δ = δ_array[1]
+δ = δ_array[2]
+δ = δ_array[3]
+δ = δ_array[4]
+δ = δ_array[5]
+δ = δ_array[6]
+δ = δ_array[7]
+δ = δ_array[8]
+δ = δ_array[9]
+δ = δ_array[10]
+δ = δ_array[11]
+
+knife_edge_constant = (γ_1/(ρ-1)) + ((r_1 - δ)/η) + max(0,(γ_1*(1-η-ρ)/(η*(ρ-1)))) - max(r_2,(γ_2+γ_1*δ_2*λ_2)/(1-δ_2))

ConstantElasticity/3-TransitionDynamics.R

@@ -0,0 +1,69 @@
+# Transition dynamics
+
+k1 <- function(t){
+  numerator = c1_forward_shooting*exp((δ-r_1)*t)
+  denom1 = ((q*w_2)/(1-q))^(ρ/(1-ρ))
+  denom2=exp((γ_1*ρ/(1-ρ))*t)
+  denomexponent = ((1-η)/ρ)-1
+  denominator = q*((q+denom1*denom2*(1-q))^denomexponent)
+  result_k1 = (numerator/denominator)^(-1/η)
+  return(result_k1)
+}
+
+k1_prime <- function(t){
+  
+  numerator = c1_forward_shooting * exp((δ-r_1)*t)
+  
+  denominator_part1 = ((q*w_2)/(1-q))^(ρ/(1-ρ))
+  denominator_part2 = exp( ( (γ_1*ρ)/(1-ρ) )*t )*(1-q)*denominator_part1
+  denominator_part3 = (q+denominator_part2)
+  denominator_part4 = denominator_part3^( ( (1-η)/ ρ ) - 1)
+  denominator_part5 = q*denominator_part4
+  denominator = denominator_part5
+
+  exponent_total = (-1/η)
+  result_k1 = (numerator/denominator)^exponent_total
+  return(result_k1)
+}
+
+
+k2 <- function(t){
+  constant1_k2 = (w_2 * λ_2 * δ_2)/(δ-r_2)
+  constant2_k2 = ((1-λ_2)/λ_2)^(δ_2*(1-λ_2))
+  constant_exponent = (1/(1-δ_2))
+  constant_term_k2 = (constant1_k2*constant2_k2)^constant_exponent
+  exponential_growth_term_k2 = (γ_2+γ_1*δ_2*λ_2)/(1-δ_2) + γ_1
+  exponential_term_k2 = exp(exponential_growth_term_k2*t)
+  result_k2 = constant_term_k2*exponential_term_k2
+  return(result_k2)
+}
+
+l1 <- function(t, k1_t, L_t){
+  constant_term_l1 = ((q*w_2)/(1-q))^(1/(ρ-1))
+  exponential_term_l1 = exp((γ_1/(ρ-1))*t)
+  result_l1 = (constant_term_l1*exponential_term_l1*k1_t)/L_t
+  return(result_l1)
+}
+
+l2 <- function(t, k2_t, L_t){
+  constant_term_l2 = (1-λ_2)/(w_2*λ_2)
+  exponential_term_l2 = exp(γ_1*t)
+  result_l2 = (constant_term_l2*k2_t) / (exponential_term_l2*L_t)
+  return(result_l2)
+}
+
+dL <- function(t,k2_t,L_t,l2_t){
+  dynamical_term_dL= ((k2_t^λ_2)*((L_t*l2_t)^(1-λ_2)))^(δ_2)
+  result_dL = r_2*L_t + β_2*exp(γ_2*t)*dynamical_term_dL
+  return(result_dL)
+}
+
+wagesPaidOrEarned <- function(L_t, t, l1_t, l2_t){
+  result_wages = L_t*w_2*exp(γ_1*t)*(1-l1_t-l2_t)
+  return(result_wages)
+}
+
+dK <- function(t,K_t,k1_t,k2_t,L_t,l1_t,l2_t){
+  result_dK = r_1*K_t - k1_t - k2_t + L_t*w_2*exp(γ_1*t)*(1-l1_t-l2_t)
+  return(result_dK)
+}

ConstantElasticity/4-ForwardShooting.R

@@ -0,0 +1,145 @@
+# Forward shooting
+options(digits=22)
+## Options
+
+## Definitions
+K_t = K_init
+L_t = L_init
+max = max(times_forward_shooting)
+l =length(times_forward_shooting)
+
+## Initialization
+pb <- txtProgressBar(min = 1, max = l, style = 3)
+K_array_forward_shooting <- vector(mode ="numeric", length=l)
+L_array_forward_shooting <- vector(mode ="numeric", length=l)
+k1_array_forward_shooting <- vector(mode ="numeric", length=l)
+k2_array_forward_shooting <- vector(mode ="numeric", length=l)
+l1_array_forward_shooting <- vector(mode ="numeric", length=l)
+l2_array_forward_shooting <- vector(mode ="numeric", length=l)
+wages_array_forward_shooting <- vector(mode ="numeric", length=l)
+
+## Calculations
+comienzo = Sys.time()
+i=1
+for(i in c(1:l)){
+  t = times_forward_shooting[i]
+
+  k1_t = k1(t)
+  k2_t = k2(t)
+  l1_t = l1(t, k1_t, L_t)
+  l2_t = l2(t, k2_t, L_t)
+  wages_t = wagesPaidOrEarned(L_t, t, l1_t, l2_t)
+  L_t = L_t + dL(t,k2_t,L_t,l2_t)*stepsize
+  K_t = K_t + dK(t,K_t,k1_t,k2_t,L_t,l1_t,l2_t)*stepsize
+  
+  ## print(c(k1_t, k2_t, l1_t, l2_t, L_t, K_t))
+  ## if(i %% 100 == 0){
+  ##  print(c(i, k2_t))
+  ##  Sys.sleep(1)
+  ## }
+  k1_array_forward_shooting[i] <- k1_t 
+  k2_array_forward_shooting[i] <- k2_t
+  wages_array_forward_shooting[i] <- wages_t
+  l1_array_forward_shooting[i] <- l1_t
+  l2_array_forward_shooting[i] <- l2_t
+  K_array_forward_shooting[i] <- K_t
+  L_array_forward_shooting[i] <- L_t
+  setTxtProgressBar(pb, i)
+  # i = i+1
+}
+fin = Sys.time()
+timetaken = fin-comienzo
+timetaken
+
+## Checking condition
+### options(digits=22) 
+### sink("/dev/null")
+
+### Shouldn't be negative
+K_array_forward_shooting[l]
+K_growth = (K_array_forward_shooting[l]-K_array_forward_shooting[l-1])/K_array_forward_shooting[l-1]/(1*stepsize) 
+K_growth
+sum(K_array_forward_shooting<0)
+
+L_array_forward_shooting[l]
+L_growth = (L_array_forward_shooting[l]-L_array_forward_shooting[l-1])/L_array_forward_shooting[l-1]/(1*stepsize)
+L_growth
+sum(L_array_forward_shooting<0)
+
+### Should be similar growth as in the model
+k1_growth = (k1_array_forward_shooting[l]-k1_array_forward_shooting[l-1])/k1_array_forward_shooting[l-1]/stepsize
+k1_growth
+
+k2_array_forward_shooting[l]
+k2_growth = (k2_array_forward_shooting[l]-k2_array_forward_shooting[l-1])/k2_array_forward_shooting[l-1]/stepsize
+k2_growth
+
+k1_array_forward_shooting[l]
+k1_array_forward_shooting[l]/K_array_forward_shooting[l]
+plot((k1_array_forward_shooting/K_array_forward_shooting)[(l-100):l])
+## sink()
+
+l2_array_forward_shooting[l/2]
+
+## Plotting (useful for back-and-forth plotting, but redundant with Plotting.R)
+K_array_plotting <- K_array_forward_shooting
+L_array_plotting <- L_array_forward_shooting
+k1_array_plotting <- k1_array_forward_shooting
+k2_array_plotting <- k2_array_forward_shooting
+l1_array_plotting <- l1_array_forward_shooting
+l2_array_plotting <- l2_array_forward_shooting
+times <- times_forward_shooting 
+
+## Knife edge study old
+
+if(knife_edge_constant>0){
+  knife_edge_constant_above = knife_edge_constant
+  K_above_knife_edge <- K_array_forward_shooting
+  L_above_knife_edge <- L_array_forward_shooting
+  k1_above_knife_edge <- k1_array_forward_shooting
+  k2_above_knife_edge <- k2_array_forward_shooting
+  wages_above_knife_edge <- wages_array_forward_shooting
+  l1_above_knife_edge <- l1_array_forward_shooting
+  l2_above_knife_edge <- l2_array_forward_shooting
+  l3_above_knife_edge <- 1- l1_array_forward_shooting - l2_array_forward_shooting
+}
+
+if(knife_edge_constant == 0){
+  knife_edge_constant_0 = knife_edge_constant ## 0
+  K_at_knife_edge <- K_array_forward_shooting
+  L_at_knife_edge <- L_array_forward_shooting
+  k1_at_knife_edge <- k1_array_forward_shooting
+  k2_at_knife_edge <- k2_array_forward_shooting
+  wages_at_knife_edge <- wages_array_forward_shooting
+  l1_at_knife_edge <- l1_array_forward_shooting
+  l2_at_knife_edge <- l2_array_forward_shooting
+  l3_at_knife_edge <- 1- l1_array_forward_shooting - l2_array_forward_shooting
+}
+
+if(knife_edge_constant <0){
+  knife_edge_constant_below = knife_edge_constant
+  K_below_knife_edge <- K_array_forward_shooting
+  L_below_knife_edge <- L_array_forward_shooting
+  k1_below_knife_edge <- k1_array_forward_shooting
+  k2_below_knife_edge <- k2_array_forward_shooting
+  wages_below_knife_edge <- wages_array_forward_shooting
+  l1_below_knife_edge <- l1_array_forward_shooting
+  l2_below_knife_edge <- l2_array_forward_shooting
+  l3_below_knife_edge <- 1- l1_array_forward_shooting - l2_array_forward_shooting
+  
+}
+
+## knife edge study
+times_array <- c(times_array,times)
+knife_edge_constant_array <- c(knife_edge_constant_array, rep(knife_edge_constant, length(times)))
+K_knife_edge <- c(K_knife_edge,K_array_forward_shooting)
+L_knife_edge <- c(L_knife_edge,L_array_forward_shooting)
+k1_knife_edge <- c(k1_knife_edge,k1_array_forward_shooting)
+k2_knife_edge <- c(k2_knife_edge,k2_array_forward_shooting)
+wages_knife_edge <- c(wages_knife_edge,wages_array_forward_shooting)
+l1_knife_edge <- c(l1_knife_edge,l1_array_forward_shooting)
+l2_knife_edge <- c(l2_knife_edge,l2_array_forward_shooting)
+l3_knife_edge <- c(l3_knife_edge,1- l1_array_forward_shooting - l2_array_forward_shooting)
+
+## For convenience
+K_array_forward_shooting[l]

ConstantElasticity/5-Plotting.R

@@ -0,0 +1,509 @@
+# Plotting
+## install.packages("tidyverse") <- Not totally necessary.
+## install.packages("ggplot2")
+
+library(ggplot2)
+library(ggsci) ## Used for scales
+
+## General variables
+saveplots=TRUE
+### saveplots=FALSE
+### shootingtype="reverse"
+
+shootingtype="forward"
+utilityfunction="ConstantElasticityOfSubstitution"
+directory = paste("/home/nuno/Documents/core/srf/fresh/simulations/rcode/", utilityfunction, "/plots/", shootingtype, "shooting", last, "years", sep="")
+
+directory
+setwd(directory)
+getwd()
+
+height = 5
+width = floor(height*(1+sqrt(5))/2)
+imagenumbercounter = 1
+saveplot = function(imagename){
+  if(saveplots){
+    ggsave(paste(imagenumbercounter,"_", imagename, "_",shootingtype, "shooting",".png", sep =""), units="in", width=width, height=height)  
+    imagenumbercounter <<- imagenumbercounter+1
+    ## https://stackoverflow.com/questions/1236620/global-variables-in-r
+  }
+}
+options(digits=1) ## Just for display
+
+## Data
+if(shootingtype=="forward"){
+  K_array_plotting <- K_array_forward_shooting
+  L_array_plotting <- L_array_forward_shooting
+  k1_array_plotting <- k1_array_forward_shooting
+  k2_array_plotting <- k2_array_forward_shooting
+  wages_array_plotting <- -wages_array_forward_shooting
+  l1_array_plotting <- l1_array_forward_shooting
+  l2_array_plotting <- l2_array_forward_shooting
+} else if(shootingtype=="reverse"){
+  K_array_plotting <- K_array_reverse_shooting
+  L_array_plotting <- L_array_reverse_shooting
+  k1_array_plotting <- k1_array_reverse_shooting
+  k2_array_plotting <- k2_array_reverse_shooting
+  l1_array_plotting <- l1_array_reverse_shooting
+  l2_array_plotting <- l2_array_reverse_shooting
+}
+
+## K and L
+
+xs <- list()
+xs$values <- c(K_array_plotting, L_array_plotting)
+xs$var <- c(rep("K", length(K_array_plotting)), rep("L", length(L_array_plotting)))
+xs$times = rep(times,2)
+xs <- as.data.frame(xs)
+title_text = "Evolution of state variables"
+(ggplot(data = xs, aes(x = times, y= values, color = var)) +
+  geom_line(size = 0.5)
+  + labs(
+    title=title_text,
+    x="Year since start", 
+    y="Capital and labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("K", "L"),
+    labels=c("Capital", "Labor")
+    
+  )
+  + scale_y_continuous(breaks = seq(min(K_array_plotting), max(K_array_plotting), length.out=5))
+)
+saveplot("StateVariablesKL")
+
+## K and L: log plot
+title_text = "Evolution of state variables"
+(ggplot(data = xs, aes(x = times, y= values, color = var)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle="(logarithmic scale)",
+    x="Year since start", 
+    y="Capital and labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5),
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("K", "L"),
+    labels=c("Capital", "Labor")
+    
+  )
+  + scale_y_continuous(trans = 'log2')
+)
+saveplot("StateVariablesKLlogplot")
+
+## Only L
+equil2 <- list()
+equil2$values <- c(L_array_plotting)
+equil2$times = times
+equil2 <- as.data.frame(equil2)
+
+title_text = "Evolution of movement size"
+(ggplot(data = equil2, aes(x = times, y= values)) 
+  +geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+)
+saveplot("StateVariableL")
+
+
+## l1, l2, s4
+sigmas <- list()
+s4_array_plotting=1-l1_array_plotting-l2_array_plotting
+sigmas$values <- c(l1_array_plotting, l2_array_plotting,s4_array_plotting)
+sigmas$fraction <- c(rep("l1", length(l1_array_plotting)), rep("l2", length(l2_array_plotting)),rep("s4", length(l2_array_plotting)))
+sigmas$times = rep(times,3)
+sigmas <- as.data.frame(sigmas)
+
+title_text = "Evolution of labor fractions"
+(ggplot(data = sigmas, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Fraction of total labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("l1", "l2", "s4"),
+    labels=c("Direct work", "Movement building", "Money-making")
+  )
+)
+saveplot("MovementFractions")
+
+## l1 and l2
+sigmas <- list()
+sigmas$values <- c(l1_array_plotting, l2_array_plotting)
+sigmas$fraction <- c(rep("l1", length(l1_array_plotting)), rep("l2", length(l2_array_plotting)))
+sigmas$times = rep(times,2)
+sigmas <- as.data.frame(sigmas)
+
+title_text = "Evolution of labor fractions"
+subtitle_text = "(only direct work and movement building)"
+(ggplot(data = sigmas, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle=subtitle_text, 
+    x="Year since start", 
+    y="Fraction of total labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("l1", "l2"),
+    labels=c("Direct work", "Movement building")
+    
+  )
+)
+saveplot("MovementFractionsl1l2")
+
+## Just direct work.
+sigma_1 <- list()
+sigma_1$values <- c(l1_array_plotting)
+sigma_1$times = times
+sigma_1 <- as.data.frame(sigma_1)
+
+title_text = "Evolution of direct work \nas a fraction of total labor"
+(ggplot(data = sigma_1, aes(x = times, y= values)) 
+  +geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Fraction of total labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+)
+saveplot("MovementFractionsl1")
+
+
+## Just movement building
+sigma_2 <- list()
+sigma_2$values <- c(l2_array_plotting)
+sigma_2$times = times
+sigma_2 <- as.data.frame(sigma_2)
+
+title_text = "Evolution of movement building\nas a fraction of total labor"
+(ggplot(data = sigma_2, aes(x = times, y= values)) 
+  +geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Fraction of labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+)
+saveplot("MovementFractionsl2")
+
+## Just direct work - logplot
+sigma_1 <- list()
+sigma_1$values <- c(l1_array_plotting)
+sigma_1$times = times
+sigma_1 <- as.data.frame(sigma_1)
+
+title_text = "Evolution of direct work \nas a fraction of labor"
+(ggplot(data = sigma_1, aes(x = times, y= values)) 
+  +geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle="(logarithmic scale)", 
+    x="Year since start", 
+    y="Fraction of labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+  + scale_y_continuous(trans = 'log2')
+)
+saveplot("MovementFractionsl1logplot")
+
+## Just movement building
+sigma_2 <- list()
+sigma_2$values <- c(l2_array_plotting)
+sigma_2$times = times
+sigma_2 <- as.data.frame(sigma_2)
+
+title_text = "Evolution of movement building\nas a fraction of labor"
+(ggplot(data = sigma_2, aes(x = times, y= values)) 
+  +geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle="(logarithmic scale)", 
+    x="Year since start", 
+    y="Fraction of total labor" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.subtitle = element_text(hjust = 0.5), 
+    plot.title = element_text(hjust = 0.5),
+    legend.position="bottom",
+    legend.box="vertical"
+  )
+  + scale_y_continuous(trans = 'log2')
+)
+saveplot("MovementFractionsl2logplot")
+
+## sigma_1*L and sigma_2*L, sigma_4*L, 
+sigmaxs <- list()
+sigmaxs$values <- c(l1_array_plotting*L_array_plotting, l2_array_plotting*L_array_plotting, s4_array_plotting*L_array_plotting)
+sigmaxs$fraction <- c(rep("l1L", length(l1_array_plotting)), rep("l2L", length(l2_array_plotting)), rep("s4L", length(l2_array_plotting)))
+sigmaxs$times = rep(times,3)
+sigmaxs <- as.data.frame(sigmaxs)
+
+title_text = "Evolution of labor\n in absolute terms"
+(ggplot(data = sigmaxs, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Absolute number\nof movement participants" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("l1L", "l2L", "s4L"),
+    labels=c("Direct work", "Movement building", "Money-making")
+    
+  )
+)
+saveplot("MovementParticipantNumbersl1l2S4")
+
+## sigma_1*L and sigma_2*L
+sigmaxs <- list()
+sigmaxs$values <- c(l1_array_plotting*L_array_plotting, l2_array_plotting*L_array_plotting)
+sigmaxs$fraction <- c(rep("l1L", length(l1_array_plotting)), rep("l2L", length(l2_array_plotting)))
+sigmaxs$times = rep(times,2)
+sigmaxs <- as.data.frame(sigmaxs)
+
+title_text = "Evolution of movement participants\nin absolute terms"
+(ggplot(data = sigmaxs, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle="(only direct workers and movement builders)", 
+    x="Year since start", 
+    y="Absolute number\nof movement participants" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("l1L", "l2L"),
+    labels=c("Direct work", "Movement building")
+    
+  )
+)
+saveplot("MovementParticipantNumbersl1l2")
+
+## sigma_1*L and sigma_2*L logplot
+sigmaxs <- list()
+sigmaxs$values <- c(l1_array_plotting*L_array_plotting, l2_array_plotting*L_array_plotting)
+sigmaxs$fraction <- c(rep("l1L", length(l1_array_plotting)), rep("l2L", length(l2_array_plotting)))
+sigmaxs$times = rep(times,2)
+sigmaxs <- as.data.frame(sigmaxs)
+
+title_text = "Evolution of movement participants\nin absolute terms \n(logarithmic scale)"
+(ggplot(data = sigmaxs, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle="(only direct workers and movement builders)", 
+    x="Year since start", 
+    y="Absolute number\nof movement participants" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("l1L", "l2L"),
+    labels=c("Direct work", "Movement building")
+    
+  )
+  + scale_y_continuous(trans = 'log2')
+)
+saveplot("MovementParticipantNumbersl1l2logplot")
+
+## alpha1 and alpha2
+alphas <- list()
+alphas$values <- c(k1_array_plotting, k2_array_plotting)
+alphas$fraction <- c(rep("k1", length(l1_array_plotting)), rep("k2", length(l2_array_plotting)))
+alphas$times = rep(times,2)
+alphas <- as.data.frame(alphas)
+
+title_text = "Evolution of spending"
+(ggplot(data = alphas, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Spending" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("k1", "k2"),
+    labels=c("Direct spending","Spending on movement building")
+    
+  )
+)
+saveplot("Spendingk1k2")
+
+## alphk1 and alphk2: log plot
+
+alphas <- list()
+alphas$values <- c(k1_array_plotting, k2_array_plotting)
+alphas$fraction <- c(rep("k1", length(l1_array_plotting)), rep("k2", length(l2_array_plotting)))
+alphas$times = rep(times,2)
+alphas <- as.data.frame(alphas)
+
+title_text = "Evolution of spending"
+(ggplot(data = alphas, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    subtitle="(logarithmic scale)",
+    x="Year since start", 
+    y="Spending" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5),
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("k1", "k2"),
+    labels=c("Direct spending","Spending on movement building")
+  )
+  + scale_y_continuous(trans = 'log2')
+)
+saveplot("Spendingk1k2logplot")
+
+## alpha1 and alpha2 and wages
+alphas <- list()
+alphas$values <- c(k1_array_plotting, k2_array_plotting, wages_array_plotting)
+alphas$fraction <- c(rep("k1", length(l1_array_plotting)), rep("k2", length(l2_array_plotting)), rep("wages", length(wages_array_plotting)))
+alphas$times = rep(times,3)
+alphas <- as.data.frame(alphas)
+
+title_text = "Evolution of spending"
+(ggplot(data = alphas, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Spending" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("k1", "k2", "wages"),
+    labels=c("Direct spending","Spending on movement building", "Spending on wages")
+  )
+)
+saveplot("Spendingk1k2k3")
+
+
+## alpha1 and alpha2 and wages: logplot
+alphas <- list()
+alphas$values <- c(k1_array_plotting, k2_array_plotting, wages_array_plotting)
+alphas$fraction <- c(rep("k1", length(l1_array_plotting)), rep("k2", length(l2_array_plotting)), rep("wages", length(wages_array_plotting)))
+alphas$times = rep(times,3)
+alphas <- as.data.frame(alphas)
+
+title_text = "Evolution of spending"
+(ggplot(data = alphas, aes(x = times, y= values, color = fraction)) +
+    geom_line(size = 0.5)
+  +labs(
+    title=title_text,
+    x="Year since start", 
+    y="Spending" 
+  )
+  +theme(
+    legend.title = element_blank(),
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.position="bottom",
+    legend.direction="vertical"
+  )
+  + scale_color_lancet(
+    breaks=c("k1", "k2", "wages"),
+    labels=c("Direct spending","Spending on movement building", "Spending on wages")
+  )
+  + scale_y_continuous(trans = 'log2')
+)
+saveplot("Spendingk1k2k3_logplot")

ConstantElasticity/6-Plotting-many.R

@@ -0,0 +1,260 @@
+library(ggplot2)
+library(ggsci) ## Used for scales
+
+## Directory
+directory = paste("/home/nuno/Documents/core/srf/fresh/simulations/rcode/", utilityfunction, "/plots/", sep="")
+directory
+setwd(directory)
+height = 5
+width = floor(height*(1+sqrt(5))/2)
+imagenumbercounter = 1
+saveplots=TRUE
+saveplot = function(imagename){
+  if(saveplots){
+    ggsave(paste(imagenumbercounter,"_", imagename, "_knifeedge",".png", sep =""), units="in", width=width, height=height)  
+    imagenumbercounter <<- imagenumbercounter+1
+    ## https://stackoverflow.com/questions/1236620/global-variables-in-r
+  }
+}
+options(digits=2) ## Just for display
+
+## Capital
+data = list()
+data$capital = c(K_above_knife_edge,K_at_knife_edge,K_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+options(digits=2) ## Just for display
+variable = "K"
+meaning = "Capital"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=capital)) +
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## Labor
+data = list()
+data$labor = c(L_above_knife_edge,L_at_knife_edge,L_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "L"
+meaning = "Labor"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=labor))+
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## l1
+data = list()
+data$l1 = c(l1_above_knife_edge,l1_at_knife_edge,l1_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "l1"
+meaning = "Direct work"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=l1))+
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## l2
+data = list()
+data$l2 = c(l2_above_knife_edge,l2_at_knife_edge,l2_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "l2"
+meaning = "Movement building labor"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=l2))+
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## l3
+data = list()
+data$l3 = c(l3_above_knife_edge,l3_at_knife_edge,l3_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "l3"
+meaning = "Money-making/hired labor"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=l3))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot("Money-makingORhired labor")
+
+## k1
+data = list()
+data$k1 = c(k1_above_knife_edge,k1_at_knife_edge,k1_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "k1"
+meaning = "Direct spending"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=k1))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+
+ggplot(data = data, aes(x=times, y=k1))+
+  geom_point(aes(color=category), size=0.1)+
+  scale_y_continuous(trans = 'log2')+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## k2
+data = list()
+data$k2 = c(k2_above_knife_edge,k2_at_knife_edge,k2_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "k2"
+meaning = "Spending on movement building"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=k2))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+ggplot(data = data, aes(x=times, y=k2))+
+  geom_point(aes(color=category), size=0.1)+
+  scale_y_continuous(trans = 'log2')+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+
+## wages
+data = list()
+data$wages = c(wages_above_knife_edge,wages_at_knife_edge,wages_below_knife_edge)
+data$category = c(rep(knife_edge_constant_above, length(times)), rep(knife_edge_constant_0, length(times)), rep(knife_edge_constant_below, length(times)))
+data$times = rep(times, 3)
+data <- as.data.frame(data)
+
+variable = "L(t)*w2*exp(γ1*t)*(1-l1-l2)"
+meaning = "wages"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=wages))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)

ConstantElasticity/7-Plotting-many.R

@@ -0,0 +1,279 @@
+library(ggplot2)
+library(ggsci) ## Used for scales
+
+## Directory
+directory = paste("/home/nuno/Documents/core/srf/fresh/simulations/rcode/", utilityfunction, "/plots/", sep="")
+directory
+setwd(directory)
+
+## Initializing elements
+times_array <- c() ## times_array <- c(rep(times, length(K_knife_edge)/length(times)))
+knife_edge_constant_array <- c() 
+K_knife_edge <- c()
+L_knife_edge <- c()
+k1_knife_edge <- c()
+k2_knife_edge <- c()
+wages_knife_edge <- c()
+l1_knife_edge <- c()
+l2_knife_edge <- c()
+l3_knife_edge <- c()
+
+## Saving plots
+height = 5
+width = floor(height*(1+sqrt(5))/2)
+imagenumbercounter = 1
+saveplots=TRUE
+saveplot = function(imagename){
+  if(saveplots){
+    ggsave(paste(imagenumbercounter,"_", imagename, "_knifeedge",".png", sep =""), units="in", width=width, height=height)  
+    imagenumbercounter <<- imagenumbercounter+1
+    ## https://stackoverflow.com/questions/1236620/global-variables-in-r
+  }
+}
+options(digits=2) ## Just for display
+
+## Capital
+data = list()
+data$capital = K_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+options(digits=2) ## Just for display
+variable = "K"
+meaning = "Capital"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=capital)) +
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  ) # +
+ #scale_y_continuous(trans = 'log2')
+
+saveplot(meaning)
+
+## Labor
+data = list()
+data$labor = L_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "L"
+meaning = "Labor"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=labor))+
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+
+saveplot(meaning)
+
+## l1
+data = list()
+data$l1 = l1_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "l1"
+meaning = "Direct work"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=l1))+
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## l2
+data = list()
+data$l2 = l2_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "l2"
+meaning = "Movement building labor"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=l2))+
+  geom_point(aes(color=category), size=0.05) +
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+## l3
+data = list()
+data$l3 = l3_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "l3"
+meaning = "Money-making/hired labor"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=l3))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot("Money-makingORhired labor")
+
+## k1
+data = list()
+data$k1 = k1_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "k1"
+meaning = "Direct spending"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=k1))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+ggplot(data = data, aes(x=times, y=k1))+
+  geom_point(aes(color=category), size=0.1)+
+  scale_y_continuous(trans = 'log2')+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(paste(meaning,"logplot", sep="_"))
+
+## k2
+data = list()
+data$k2 = k2_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "k2"
+meaning = "Spending on movement building"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=k2))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)
+
+ggplot(data = data, aes(x=times, y=k2))+
+  geom_point(aes(color=category), size=0.1)+
+  scale_y_continuous(trans = 'log2')+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(paste(meaning,"logplot", sep="_"))
+
+## wages
+data = list()
+data$wages = wages_knife_edge
+data$category = knife_edge_constant_array
+data$times = times_array
+data <- as.data.frame(data)
+
+variable = "L(t)*w2*exp(γ1*t)*(1-l1-l2)"
+meaning = "wages"
+title_text = paste("Evolution of ", tolower(meaning))
+ggplot(data = data, aes(x=times, y=wages))+
+  geom_point(aes(color=category), size=0.1)+
+  labs(
+    title=title_text,
+    subtitle=paste("(",variable,")",sep=""),
+    color="Knife edge",
+    x="Year since start", 
+    y=meaning
+  ) +
+  theme(
+    plot.title = element_text(hjust = 0.5),
+    plot.subtitle = element_text(hjust = 0.5), 
+    legend.direction="vertical"
+  )
+saveplot(meaning)