######################################################################################################################################################## # This program describes the dynamics of a tree nursery selling trees that undergo 4 growth stages before commercialisation. The trees grow from seeds # # to medium trees. Then it calculates the value of several economic indicators dependent on demand variability over time. # ######################################################################################################################################################## Year <- numeric() TotalYears <- 44 #Total number of years the simulations runs over # create arrays called for economic indicators seed<-array(NA,c(TotalYears+1,4)) # number of seeds planted each year demand <- array(NA,c(TotalYears+1,4)) # demand of trees of each stage every year population<-array(NA,c(TotalYears+1,4)) # age structured tree grown at each year where N(t+1)=L*N(t) production<-array(NA, c(TotalYears+1,4)) # trees produced minus sales per year pop3Sdemand<-array(NA, c(TotalYears+1,4)) # demand at stage 3 pop4Sdemand<-array(NA, c(TotalYears+1,4)) # demand at stage 4 imports <- array(NA,c(TotalYears+1,4)) # imports of trees of each stage per year cost <- array(NA,c(TotalYears+1,4)) # total costs per year profit <- array(NA,c(TotalYears+1,4)) # total gross margin per year diseaseProb <- array(NA,c(TotalYears+1,4)) # probability of importing a disease mean.demandS3 <- 2860 mean.demandS4 <- 715 seed.pop <- 3600 #number of seed planted each year p<-0.0001 # we assume one in 10000 of foreign trees has a disease # demand.var <-50 #demand variability # Forestry commission newly planted oak trees from 1976-2017 (estimate in thousand trees) demand.varS3 <- 0.8*c(0,0,0,450,450,390,480,600,600,570,690,750,810,1050,1440,1980,1950,3090,3510,3450,4410,5010,4320,3720,4020,3870,3960,4380,5070,3870,3720,3660,3870,3330,3690,3210,2790,2550,3300,3780,3690,4470,4410,2190,2010) demand.varS4 <- 0.2*c(0,0,0,450,450,390,480,600,600,570,690,750,810,1050,1440,1980,1950,3090,3510,3450,4410,5010,4320,3720,4020,3870,3960,4380,5070,3870,3720,3660,3870,3330,3690,3210,2790,2550,3300,3780,3690,4470,4410,2190,2010) # Initial values seed[1,] <- c(seed.pop,0,0,0) demand[1,] <- c(0,0,0,0) population[1,] <- c(seed.pop,0,0,0) pop3Sdemand[1,] <- c(0,0,0,0) pop4Sdemand[1,] <- c(0,0,0,0) imports[1,] <- pmax(0,demand[1,]-population[1,]) production[1,] <- c(seed.pop,0,0,0) cost[1,] <- rbind(0.035,0,0,0)*(seed[1,])+rbind(0.0,0.0,0.045,0.0)*(production[1,]) + rbind(0,0,0.12,0.15)*(pmax(0,imports[1,]))+ rbind(0,0,0.05,0.1)*(pmax(0,thrown[1,])) pop_sale[1, ]<-c(0,0,0,0) profit[1,] <- rbind(0,0,0.2,0.25)*demand[1,] - cost[1,] diseaseProb[1,] <- (1-((1-p)^{imports[1,]})) qij<-1 #percentage of trees moving to the next stage*survival percentage qii<-0 #percentage of trees staying in the same stage q44<-0 #determines percentage kept after a sale Lefkovitch <- matrix(c(qii,0,0,0, 0.9,qii,0,0, 0,0.8,qii,0, 0,0,0.25,q44), nrow = 4, ncol = 4, byrow = TRUE) #0.25 of 0.8 is equivalent to 0.2 of 1 for(Year in 1:TotalYears){ seed[Year+1,]<- c(seed.pop,0,0,0) #number of seed planted each year is constant #variability in the demand is introduced through a uniform or normal distribution with random variability in each stage (can be modified to other distributions) #demand[Year+1,] <- c(0,0,0,(max(0,runif(1,mean.demandS4-demand.var,mean.demandS4+demand.var))),0) #uniform distribution #demand[Year+1,] <- c(0,0,0,(max(0,rnorm(1,mean.demandS4,(0.75*demand.var))))) #normal distribution pop3Sdemand[Year+1,]<-c(0,0,demand.varS3[Year+1],0) #demand from stage 3 pop4Sdemand[Year+1,]<-c(0,0,0,demand.varS4[Year+1]) #demand from stage 4 demand[Year+1,] <-c(0,0,pop3Sdemand[Year+1,3],pop4Sdemand[Year+1,4]) population[Year+1,] <- pmax(0,(Lefkovitch)%*%(population[Year,]) + seed[Year+1,]) production[Year+1,]<-pmax(0,population[Year+1,]-demand[Year+1,]) imports[Year+1,]<-pmax(0, demand[Year+1,]-population[Year+1,]) diseaseProb[Year+1,] <- 1-((1-p)^{imports[Year+1,]}) cost[Year+1,] <- rbind(0.035,0,0,0)*(seed[Year+1,])+rbind(0.0,0.0,0.065,0.0)*(population[Year+1,]) + rbind(0,0,0.15,0.2)*(pmax(0,imports[Year+1,]))#+ rbind(0,0,0.04,0.075)*(pmax(0,thrown[Year,])) profit[Year+1,] = rbind(0,0,0.4,0.45)*demand[Year+1,] - cost[Year+1,] }