###################################################################################################################### # This program describes the dynamics of an economic model for a tree nursery with trees moving along 4 stages, from # # seeds to medium trees. It calculates the average gross margin and the probability of introducing a forest disease # # dependent on planting rate and demand variability of trees in stage 4, assuming sales in stage 4 only. This can # # be modified to include Stage 3. Here we use data of newly planted oak trees from the Forestry Commission as an # # indicator for tree planting rates and approximate costs of produced and imported saplings from nursery growers. # # ###################################################################################################################### ### Initial conditions ### Year <- numeric() TotalYears <- 45 # Total number of years the simulation runs over n <- numeric() TotalRuns <- 100 # Number of runs to obtain average dynamics M<- numeric() Mseq = seq(400,1500,10) # Range of planting rates (thousand trees) d <- numeric() demand.var <- seq(0,500,5) # Range of demand variability (thousand trees) k <- numeric() TotalK <- 0.15 #production costs p <- 0.0001 Stages <- 4 ### create arrays for system dynamics ### ### Gross margin ### n_ProfitMean <- array(NA,c(TotalYears+1,Stages)) # Gross margin mean of each seed planting rate for all the runs i.e. dataset of TotalYears*Stages n_StageMean <- array(NA,c(length(Mseq),Stages)) # Gross margin mean of each stage after n runs i.e. vector of length(Stages) M_StageMean <- array(NA,c(length(Mseq),Stages)) # Gross margin mean of each stage after n runs for all planting rates i.e. a dataset of length(Mseq)*Stages ### Costs ### n_CostMean <- array(NA,c(TotalYears+1,Stages)) n_StageCostMean <- array(NA,c(length(Mseq),Stages)) M_StageCostMean <- array(NA,c(length(Mseq),Stages)) ### Imports ### n_ImportsMean <- array(NA,c(TotalYears+1,Stages)) n_StageImportsMean <- array(NA,c(length(Mseq),Stages)) M_StageImportsMean <- array(NA,c(length(Mseq),Stages)) ### Probabilities of importing a disease ### n_Prob <- array(NA,c(TotalYears+1,Stages)) n_StageProbMean <- array(NA,c(length(Mseq),Stages)) M_Prob<- array(NA,c(length(Mseq),Stages)) ### Final calculated values ### AllStageProfit <- array(NA,c(length(Mseq),Stages)) # Gross margin sum of means of all stages after n runs for each planting rate and fixed demand variability AllStageCost <- array(NA,c(length(Mseq),Stages)) AllStageImports <- array(NA,c(length(Mseq),Stages)) AllStageProbs <- array(NA,c(length(Mseq),Stages)) ### Create lists of variables for each planting rate and demand variability value ### Profit_var <- numeric() Cost_var <- numeric() Imports_var <- numeric() Prob_var <- numeric() ### create arrays of tree grownt and economic factors ### seed <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # number of seeds planted for a set of runs over time. A vector of (N,0,...,0) demand <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # demand of trees of each class every year. population <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # population is the age structured population at each year so that N(t+1)=L*N(t) production<-array(NA,c(TotalYears+1,Stages,TotalRuns)) # production is the age structured population at each year undergoing management #pop3Sdemand <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # variable demand of trees in stage3 every year pop4Sdemand <- array(NA, c(TotalYears+1,Stages,TotalRuns)) # variable demand of trees in stage4 every year imports <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # imports of trees of each class every year cost <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # total costs (production, imports and stock holding) over time profit <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # total profit over time pop_sale <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # trees left after a sale diseaseProb <- array(NA,c(TotalYears+1,Stages,TotalRuns)) # probability of importing a disease max_profit <- array(NA,c(length(demand.var),Stages)) #maximum profit for each demand variability value over all planting rates opt_planting <- array(NA,c(length(demand.var),Stages)) #to which planting rate the maximum profit belongs ### Leftkovitch tree growth model: age structured matrix ### T_BS <- 1.0 T_SL <- 1.0 T_LM <- 1.0 T_MM <- 1.0 A_B <- 1.0 A_S <- 1.0 A_L <- 1.0 A_M <- 1.0 Transition <- matrix(c(0,0,0,0, T_BS,0,0,0, 0,T_SL,0,0, 0,0,T_LM,0), nrow = 4, ncol = 4, byrow = TRUE) NoTransition <- matrix(c(1-T_BS,0,0,0, 0,1-T_SL,0,0, 0,0,1-T_LM,0, 0,0,0,1-T_MM), nrow = 4, ncol = 4, byrow = TRUE) Survival <-matrix(c(A_B,0,0,0, 0,A_S,0,0, 0,0,A_L,0, 0,0,0,A_M),nrow = 4, ncol = 4, byrow = TRUE) Lefkovitch <-(Transition+NoTransition)%*%Survival Lefk_Trans <-(Transition)%*%Survival Lefk_NoTrans <-(NoTransition)%*%Survival ### mean demands ### #mean.demandS3 <- 4200 # mean demand stage4/year mean.demandS4 <- 1000 # mean demand stage4/year seed.pop <- 1000 # number of planted seeds/year ################################################## Running ################################################ #Runs for demand variability sequence for(d in 1:length(demand.var)){ #Runs for different planting rates for(M in 1:length(Mseq)){ #Runs economic variables n times to obtain a mean economics dynamics for(n in 1:TotalRuns){ #Runs economic variables over a total number years for(Year in 1:TotalYears){ ### initial values seed[1,,n] <- Mseq[M]*c(1,0,0,0) demand[1,,n] <- c(0,0,0,mean.demandS4) population[1,,n] <- Mseq[M]*c(1,0,0,0) pop_sale[1,,n] <- c(0,0,0,0) production[1,,n]<-Mseq[M]*c(1,0,0,0) pop4Sdemand[1,,n] <- c(0,0,0,mean.demandS4) #pop3Sdemand[1,,n] <- c(0,0,mean.demandS3,0) imports[1,,n] <- c(0,0,0,mean.demandS4) cost[1,,n] <- 0.5*(rbind(0.015,0,0,0)*(seed[1,,n])+rbind(0.005,0.02,0.035,0.0)*(production[1,,n])) + rbind(0,0,0,0.15)*(pmax(0,imports[1,,n])) profit[1,,n] <-rbind(0,0,0,0.25)*demand[1,,n] - cost[1,,n] ### number of seed planted over time is constant for each set of runs. It increases with the planting rate M after each set of runs seed[Year+1,,n] <- Mseq[M]*c(1,0,0,0) ### variability in the demand is introduced through a uniform distribution (or normal) with random variability in the demand stage ### ### Sales can happen in stages 3 and 4 or only in stage 4 ### #pop3Sdemand[Year+1,,n] <- c(0,0,0,(max(0, runif(1,mean.demandS3-demand.var[d],mean.demandS3+demand.var[d])))) #pop4Sdemand[Year+1,,n] <- c(0,0,0,(max(0, rnorm(1,mean.demandS4,(0.75*demand.var[d]))))) pop4Sdemand[Year+1,,n] <- c(0,0,0,(max(0, runif(1,mean.demandS4-demand.var[d],mean.demandS4+demand.var[d])))) #pop4Sdemand[Year+1,,n] <- c(0,0,0,(max(0, rnorm(1,mean.demandS4,(0.75*demand.var[d]))))) #demand[Year+1,,n] <- c(0,0,pop3Sdemand[Year+1,3,n],pop4Sdemand[Year+1,4,n]) demand[Year+1,,n] <- c(0,0,0,pop4Sdemand[Year+1,4,n]) ### Population increases each year in the different stages by planting seeds and transitioning stages. ### It changes according to the survival of trees in each stage and the probability of trees in each stage moving to the next stage population[Year+1,,n] <- pmax(0,(Lefkovitch)%*%(population[Year,,n]) + seed[Year+1,,n]) #population produced before sales pop_sale[Year+1,,n] <- pmax(0,pmax(0,(Lefkovitch)%*%(pop_sale[Year,,n])) + seed[Year+1,,n]-demand[Year+1,,n]) #population left after sale ### Production are the plants produced over time in each stage production[Year+1,,n] <- pmax(0,pmax(0,(Lefk_Trans)%*%production[Year,,n]+Lefk_NoTrans%*%pop_sale[Year,,n]+seed[Year+1,,n])-(Lefk_Trans)%*%demand[Year,,n]) imports[Year+1,,n] <- pmax(0, demand[Year+1,,n]-production[Year+1,,n]) cost[Year+1,,n] <- 0.5*(rbind(0.015,0,0,0)*(seed[Year+1,,n])+rbind(0.005,0.02,0.035,0.0)*(production[Year+1,,n])) + rbind(0,0,0,0.15)*(pmax(0,imports[Year+1,,n])) profit[Year+1,,n] <- rbind(0,0,0,0.25)*demand[Year+1,,n] - cost[Year+1,,n] diseaseProb[Year+1,,n] <- 1-((1-p)^{imports[Year+1,,n]}) } n_ProfitMean <- apply(profit,c(1,2),mean) n_StageMean <- colMeans(n_ProfitMean) n_CostMean <- apply(cost,c(1,2),mean) n_StageCostMean <- colMeans(n_CostMean) n_ImportsMean <- apply(imports[5:TotalYears,,],c(1,2),mean) n_StageImportsMean <- colMeans(n_ImportsMean) n_Prob <- apply(diseaseProb[5:TotalYears,,],c(1,2),mean) n_StageProbMean <- colMeans(n_Prob) } M_StageMean[M,] <- n_StageMean M_StageCostMean[M, ] <- n_StageCostMean M_StageImportsMean[M,] <- n_StageImportsMean M_Prob[M,] <- n_StageProbMean } AllStageProfit = rowSums(M_StageMean) Profit_var = c(Profit_var, AllStageProfit) AllStageCost = rowSums(M_StageCostMean) Cost_var = c(Cost_var, AllStageCost) AllStageImports =rowSums(M_StageImportsMean) Imports_var = c(Imports_var,AllStageImports) AllStageProbs = rowSums(M_Prob) Prob_var = c(Prob_var,AllStageProbs) max_profit[d,]<-apply(M_StageMean, MARGIN=c(2), max) opt_planting[d,]<-apply(M_StageMean,MARGIN=c(2),which.max) } ############################ Differential of profit respect to planting rate ######################### PROFIT <- matrix(Profit_var,ncol=length(demand.var),nrow=length(Mseq),dimnames=list(Mseq,demand.var)) row.names<-demand.var col.names<-Mseq COST <- matrix(Cost_var,ncol=length(demand.var),nrow=length(Mseq),dimnames=list(Mseq,demand.var)) row.names<-demand.var col.names<-Mseq PROFIT_df <- data.frame(x = rep(seq_len(ncol(PROFIT)), each = nrow(PROFIT)), y = rep(seq_len(nrow(PROFIT)), times = ncol(PROFIT)), z = c(PROFIT)) COST_df <- data.frame(x = rep(seq_len(ncol(PROFIT)), each = nrow(PROFIT)), y = rep(seq_len(nrow(PROFIT)), times = ncol(PROFIT)), z = c(COST)) #smooth function out so we can take the differential of profit with respect to planting rate require("mgcv") #mod <- gam(z ~ te(x, y), data = PROFIT_df) mod <- gam(z ~ te(x, y, k=c(25,50)), data = PROFIT_df) m2 <- matrix(fitted(mod), ncol = length(demand.var)) require("lattice") wireframe(m2) a<-apply(m2, 2, diff)