# Author: Ulises Diéguez-Aranda (2016)
# Examples of Clutter et al. (1983, tables 8.1 & 8.2, figures 8.1 & 8.2, pp. 212-228)

# The growth of even-aged timber stands - volumetric considerations ------------
EstimateVolume <- function(t){  # Per-acre yiedl un cunits at the end of the rotation 
        Yt <- 100 * (1 - exp(-0.05 * t)) ^ 2
        return(Yt)
}
age <- 10:40
yield <- EstimateVolume(age)
mai <- yield / age
cai <- c(yield[-1] - yield[-length(yield)], NA)
table.8.1 <- data.frame(age, yield, mai, cai)
write.csv(table.8.1, file = "table_8_1_ClutterEtlAl1983.csv", row.names = FALSE)

# The economic growth of even-aged timber stands - economic considerations -----
CalculateNFV <- function(S, Vt, R, P, i, A, t){  # Per-acre net return at harvest
  # S = Per-cunit stumpage price
  # Vt = Per-acre yield in cunits at the end of the rotation
  # R = Rotation length
  # P = Per-acre regeneration cost
  # i = Inflation-free interest rate
  # A = Annual per-acre ad valorem tax and administration cost
  NFV <- S * Vt - P * (1 + i) ^ R - A * ((1 + i) ^ R - 1) / i
  return(NFV)
}
CalculateNPV <- function(NFV, R, i){  # Per-acre net present value
  # R = Rotation length
  # i = Inflation-free interest rate
  NPV <- NFV / (1 + i) ^ R
  return(NPV)
}
CalculateLEV <- function(NFV, R, i){  # Per-acre land expectation value
  # R = Rotation length
  # i = Inflation-free interest rate
  LEV <- NFV / ((1 + i) ^ R - 1)
  return(LEV)
}
nfv <- CalculateNFV(S = 35, Vt = yield, R = age, P = 100, i = 0.04, A = 1.5)
npv <- CalculateNPV(NFV = nfv, R = age, i = 0.04)
lev <- CalculateLEV(NFV = nfv, R = age, i = 0.04)
table.8.2 <- data.frame(age, yield, mai, cai, nfv, npv, lev)
write.csv(table.8.2, file = "table_8_2_ClutterEtlAl1983.csv", row.names = FALSE)
summary(mai)
summary(cai)

# Yield graph
pdf(file = "Graph_Yield_UDA.pdf", width = 5, height = 6)
par(mar = c(5, 5, 2, 5))
plot(age, yield, type = "l", xlab = "Age (yr)", ylab = "Yield (cunits/acre)")
dev.off()

# Yield, MAI, and CAI relationship
pdf(file = "Graph_Yield_MAI_CAI_UDA.pdf", width = 5, height = 6)
par(mar = c(5, 5, 2, 5))
plot(age, yield, type = "l", xlab = "Age (yr)", ylab = "Yield (cunits/acre)")
par(new = T)
plot(age, mai, type = "l", lty = 2, axes = F, xlab = NA, ylab = NA, ylim = c(1, 2.5))
par(new = T)
plot(age, cai, type = "l", lty = 3, axes = F, xlab = NA, ylab = NA, ylim = c(1, 2.5))
axis(side = 4)
mtext(side = 4, line = 3, 'Growth (cunits/acre/year)')
legend("bottom", bty = "n", legend = c("Yield", "MAI", "CAI"), lty = c(1, 2, 3))
dev.off()

# Yield and LEV relationship
pdf(file = "Graph_Yield_LEV_UDA.pdf", width = 5, height = 6)
par(mar = c(5, 5, 2, 5))
plot(age, yield, type = "l", xlab = "Age (yr)", ylab = "Yield (cunits/acre)")
par(new = T)
plot(age, lev, type = "l", lty = 2, axes = F, xlab = NA, ylab = NA)
axis(side = 4)
mtext(side = 4, line = 3, 'LEV ($/acre)')
legend("bottom", bty = "n", legend = c("Yield", "LEV"), lty = c(1, 2))
dev.off()

# Decisions concerning existing stands -----------------------------------------
years.hence <- 0:10
volume.hence <- c(10.1, 11.8, 13.7, 15.2, 16.7, 18.1, 19.3, 20.3, 21.1, 21.8, 22.5)
table.8.5 <- data.frame(years.hence, volume.hence)

CalculateForestValue <- function(S, Vt, n, i, MaxLEV, A){  # Per-acre net return at harvest
  # S = Per-cunit stumpage price
  # Vt = Per-acre initial harvest volume in cunits
  # n = Years until initial harvest
  # i = Inflation-free interest rate
  # A = Annual per-acre ad valorem tax and administration cost
  PVIH <- S * Vt / (1 + i) ^ n  # Present value of initial harvest
  PVLEV <- MaxLEV  / (1 + i) ^ n  # Present value of LEV
  PVTAC <- - A * ((1 + i) ^ n - 1) / (i * (1 + i) ^ n)  # Present value of taxes and administrative costs
  ForestValue <- PVIH + PVLEV + PVTAC
  return(data.frame(n, Vt, PVIH, PVLEV, PVTAC, ForestValue))
}

max.lev <- max(lev)  # Clutter et al. (1983) used a value of 967.73
table.8.7 <- CalculateForestValue(S = 35, Vt = volume.hence, n = years.hence, i = 0.04, MaxLEV = max.lev, A = 1.5)
write.csv(table.8.7, file = "table_8_7_ClutterEtlAl1983.csv", row.names = FALSE)
OptimalDecision <- table.8.7[which.max(table.8.7$ForestValue), ]
TimberValue <- table.8.7$ForestValue - max.lev
IncrementalGrowingStock <- TimberValue - TimberValue[1]