library(truncnorm)
library(ggplot2)
n <- 100000
t <- 100
a <- 50
unprofit2 <- matrix(nrow = a, ncol = t)
cnt <- matrix(ncol = t)
##n is the number of firms
while (t > 0) {
pm$"profit" <- matrix(rgamma(n, shape = 2.78, rate = 20.29))
##creates gamma distribution for profit rates
Inorm <- rtruncnorm(n, a = 0, mean = 11.2, sd =  12.4)
##creates truncated normal distribution of capital intensity
pm$"CI" <- (.454/pm$"profit") * Inorm
##derives capital instensity from profit rate and the distribution function
##see paper "Labour Values, Prices of Production...Germany"
## https://www.tu-chemnitz.de/wirtschaft/vwl2/downloads/paper/froehlich/deviation.pdf
pm$"Ptotal" <- (pm$"CI"+1)*pm$"profit"
##multiply capital intensity times profit to get total profit where w = 1
rwg <- c(1:a)
rwg <- rwg/1000
rn <- a
##real wage growth test cases
unprofit <- matrix(nrow = rn)
unprofit1 <- matrix(ncol = rn, nrow = n)
unprofit[1:rn,] <- 0
##creates matrix for unprofitable firms
while (rn > 0)
  ##loop for every test case
{unprofit1[,rn] <- pm$"Ptotal" - rwg[rn];
##fills matrix column with new profit rates after a given raise in real wages
for (val in unprofit1[,rn]) {
  if (val <= 0) {unprofit[rn,] <- unprofit[rn,] + 1}};
##checks how many profit rates are now zero or lower
rn <- rn - 1;}
##moves to next test case
unprofit2[,t] <- unprofit
pm$"total" <- pm$"Ptotal" + pm$"CI"
pm$"pt1" <- pm$"Ptotal"
totals <- sum(pm$"total") * .0135
random <- c(1:500)
sample(1:n, 1)
while (totals > 0) {
  s <- sample(1:n, 1);
  r <- sample(random, 1)/1000;
  d <- pm$"pt1"[s]
  pm$"pt1"[s] <- pm$"pt1"[s] - r;
  if (pm$"pt1"[s] < 0 ) {totals <- totals - d;
  pm$"pt1"[s] <- 0
  } else {totals <- totals - d}
}
cn <- 0
for (v in pm$"pt1")
{if (v <= 0) {cn <- cn + 1}}
cnt[,t] <- cn
print(t)
t <- t - 1
}
j <- a
while (j > 0) {
unprofit[j,] <- mean(unprofit2[j,]);
j <- j - 1
}
unprofit <- (unprofit/n) * 100
unprofit[35,]
unprofit[40,]
mean(cnt)/n
u35 <- matrix(nrow = a)
u35[1:a,] <- unprofit[35,]
u40 <- matrix(nrow = a)
u40[1:a,] <- unprofit[40,]
ucn <- matrix(nrow = a)
ucn[1:a,] <- mean(cnt)/n
ucn <- ucn +u35
##turns number unprofitable into percent of firms
plot(rwg, unprofit, type = "l", xlab = "wage growth", ylab = "% firms unprofitable"
      )
lines(rwg, u35, type = "l", col = "blue")
lines(rwg, u40, type = "l", col = "red")
lines(rwg, ucn, type = "l", col = "green")
legend(.001, .0095, legend=c("correlation", "critical threshold", "current level", "current level + trade war shock"),
       col=c("black", "red", "blue", "green"), lty=1:2, cex=0.8)