setwd("~/Dropbox/Tillage x Fertility trial/Economics paper/Data/final_data_files/") ## user to define file structure yld <- read.csv("Master_Yield_1970-2014_econ.csv") yld <- yld[,1:15] ### cut unnecessary columns yld$yield_bu_ac <- as.character(yld$yield_bu_ac) yld$yield_bu_ac <- as.numeric(yld$yield_bu_ac) yld$yield_kg_ha <- as.character(yld$yield_kg_ha) yld$yield_kg_ha <- as.numeric(yld$yield_kg_ha) ### crop prices, $/bu, taken from USDA NASS quick stats corn_pr <- read.csv("corn_price_USDANASS.csv") corn_pr <- corn_pr[which(corn_pr$Period=="MARKETING YEAR"),] corn_pr <- corn_pr[,c("Year", "Value")] colnames(corn_pr)[2] <- c("Corn_pr_nominal") soy_pr <- read.csv("soy_price_USDANASS.csv") soy_pr <- soy_pr[,c("Year", "Value")] colnames(soy_pr)[2] <- c("Soy_pr_nominal") prices <- merge(corn_pr, soy_pr, by="Year", all=TRUE) #### reconstruct PPITW to update prices ppitw1990 <- read.csv("PPITW_2013-1975-1990index.csv") ppitw1990 <- ppitw1990[, c("Year", "Value")] colnames(ppitw1990)[2] <- c("PPITW90") ppitw2011 <- read.csv("PPITW_2014-1990-2011index.csv") ppitw2011 <- ppitw2011[, c("Year", "Value")] colnames(ppitw2011)[2] <- c("PPITW11") ### 1967-indexed PPITW to get index to reach back from 1970-1975 (values copied from scanned pdfs of NASS annual price summaries, 1971-1980) ppitw1967 <- data.frame(c(114, #1970 120, 127, 145, 168, 185, 192, 202, 219, 250, 280)) #1980 ppitw1967$Year <- seq(1970,1980, by=1) colnames(ppitw1967)[1] <- c("PPITW67") ppitw_comp <- merge(ppitw1967, ppitw1990, by="Year", all = TRUE) ppitw_comp <- merge(ppitw_comp, ppitw2011, by="Year", all=TRUE) ppitw_comp$PPITWcomp <- ppitw_comp$PPITW11 ### back-calculate 2011-indexed PPITW to from 1989 to 1975 using the 1990-indexed PPITW ppitw_comp$PPITWcomp[which(ppitw_comp$Year<=1989 & ppitw_comp$Year>=1975)] <- 100*(ppitw_comp$PPITW90[which(ppitw_comp$Year<=1989 & ppitw_comp$Year>=1975)]/ppitw_comp$PPITW90[44]) ### back-calculate 2011-indexed PPITW from 1974 to 1970 using the 1967-indexed PPITW (use 1980 PPITWcomp as fulcrum) ppitw_comp$PPITWcomp[which(ppitw_comp$Year<=1974 & ppitw_comp$Year>=1970)] <- ppitw_comp$PPITWcomp[11]*(ppitw_comp$PPITW67[which(ppitw_comp$Year<=1974 & ppitw_comp$Year>=1970)]/ppitw_comp$PPITW67[11]) prices <- merge(prices, ppitw_comp[,c("Year", "PPITWcomp")], by="Year", all=TRUE) ### Fertilizer costs collation ### Ammonium nitrate, urea, and ag. lime prices from USDA NASS amm_nit <- read.csv("Amm_nit_2014_1970_price.csv") amm_nit <- amm_nit[,c("Year", "Value")] colnames(amm_nit)[2] <- c("Amm_nit_pr_nominal") urea <- read.csv("Urea_2014_1970_price.csv") urea <- urea[,c("Year", "Value")] colnames(urea)[2] <- c("urea_pr_nominal") lime <- read.csv("Ag_lime_2014-1970_price.csv") lime <- lime[,c("Year", "Value")] colnames(lime)[2] <- c("Lime_pr_nominal") tab7 <- read.csv("Table7_Fert_prices_1960.csv") ### table 7 from a separate USDA report, provides clearer data on KCl and TSP prices colnames(tab7)[1] <- c("Year") tab7 <- tab7[,c("Year", "X.8", "X.10")] colnames(tab7)[2:3] <- c("TSP_pr_nominal", "KCl_pr_nominal") tab7[,2] <- as.character(tab7[,2]) tab7[,3] <- as.character(tab7[,3]) tab7 <- tab7[6:112,] tab7$Year <- as.character(tab7$Year) tab7$Year <- as.numeric(tab7$Year) tab7 <- tab7[is.finite(tab7$Year),] tab7[,2] <- as.numeric(tab7[,2]) tab7[,3] <- as.numeric(tab7[,3]) tab7 <- rbind(tab7, c(2014, 621, 601)) ### 2014 US total annual average values, as reported in NASS Quick Stats tab7 <- tab7[which(tab7$Year>=1970),] fert <- merge(amm_nit, urea, by="Year", all=TRUE) fert <- merge(fert, tab7, by="Year", all=TRUE) fert <- merge(fert, lime, by="Year", all=TRUE) prices <- merge(prices, fert, by="Year", all=TRUE) ### labor prices ####read in unsorted file from USDA NASS quick stats (filters: Economics/Expenses/Labor/Wage Rate/2015:1978) labor <- read.csv("Labor_Unsorted.csv") labor.nat <- labor[which(labor$Geo.Level=="NATIONAL"),] labor.nat.ann <- labor.nat[which(labor.nat$Period=="YEAR"),] pick <- labor.nat.ann$Data.Item[1] labor.nat.ann.type <- labor.nat.ann[which(labor.nat.ann$Data.Item==pick),]### Type = Labor, hired, wage rate, measured in $/hour wage_90_14 <- labor.nat.ann.type[1:25,c("Year", "Value")] ## select 2014-1990 wage_90_14$Value <- as.character(wage_90_14$Value) wage_90_14$Value <- as.numeric(wage_90_14$Value) #### manually read in average annual labor wages from old .wk1 files from USDA. wage_70_89 <- c(1.64, ### 1970, "per hour, wages only" wage from 70-73 1.73, 1.84, 2.00, 2.25, ### 1974, "Machine operator" wage from 74-79 2.43, 2.66, 2.87, 3.09, 3.39, 3.65, ### 1980, "all workers" wage from 80-89 4.02, 4.00, 4.11, 4.36, 4.42, 4.70, 4.87, 5.02, 5.36) wage_70_89 <- data.frame(cbind(wage_70_89, seq(1970,1989, by=1))) colnames(wage_70_89)[1:2] <- c("Value", "Year") wages <- c(wage_70_89$Value, wage_90_14$Value) years <- c(wage_70_89$Year, wage_90_14$Year) wages_comp <- data.frame(years, wages) wages_comp <- wages_comp[order(wages_comp$years),] colnames(wages_comp) <- c("Year", "Wage_pr_nominal") prices <- merge(prices, wages_comp, by="Year", all=TRUE) #### Diesel fuel price, 1970-2000, from scans of individual January Agricultural Prices reports ### 2001-2014 from USDA NASS quick stats (filters: Economics/Prices paid/fuels/Fuels, Diesel, price paid measured in $/gal) diesel <- data.frame(c(0.181, #1970 0.190, 0.190, 0.227, 0.365, 0.402, 0.413, 0.446, 0.466, 0.725, 0.977, 1.148, 1.115, 1.015, 0.998, 0.962, 0.728, 0.712, 0.729, 0.761, 0.944, 0.872, 0.816, 0.820, 0.767, 0.766, 0.920, 0.874, 0.740, 0.728, 1.080, 1.080, 0.964, 1.240, 1.310, 1.968, 2.275, 2.430, 3.619, 1.688, 2.540, 3.533, 3.740, 3.570, 3.537)) #2014 diesel$Year <- seq(1970,2014, by=1) colnames(diesel)[1] <- c("Diesel_pr_nominal") prices <- merge(prices, diesel, by="Year", all=TRUE) ### Herbicide Price reconstruction ### Glyphosate price ### Herbicide prices paid index (pegged to 1977) as read from scans of Agricultural Prices reports herb_pr <- matrix(9999, ncol=2, nrow=length(unique(yld$Year2))) herb_pr[1:45,1] <- seq(1970, 2014, by=1) herb_pr[1:22,2] <- c(62, ## 1970 63, 65, 67, 76, 102, 111, 100, 94, 96, 102, 111, 119, 125, 128, 128, 127, 124, 127, 132, 139, 151) ## 1991 herb_pr <- data.frame(herb_pr) colnames(herb_pr)[1:2] <- c("Year", "Herb_pr_ind_1977") ### read in prices taken from USDA NASS quick stats gly_14_01 <- read.csv("Gly_2014_2001_prices.csv") gly_14_01 <- gly_14_01[,c("Year", "Value")] herb_pr <- merge(herb_pr, gly_14_01, by="Year", all.x=TRUE) colnames(herb_pr)[3] <- c("Glyph4_per_gal") ### Prices found for "Glyphosate, #4 cost per gal" read from Agricultural prices scanned pdfs 91-2000 herb_pr[which(herb_pr$Year>=1991 & herb_pr$Year<=2000), c("Glyph4_per_gal")] <- c(55.40, ##1991 44.00, 52.10, 53.60, 54.10, 55.70, 56.70, 56.30, 45.50, 43.30) ##2000 ### Composite Glyphosate price history using Herbicide Price Index, 1991 price as fulcrum herb_pr$Gly_pr_per_gal_nominal <- herb_pr$Glyph4_per_gal herb_pr$Gly_pr_per_gal_nominal[which(herb_pr$Year>=1970 & herb_pr$Year<=1990)] <- (herb_pr$Herb_pr_ind_1977[which(herb_pr$Year>=1970 & herb_pr$Year<=1990)]/herb_pr$Herb_pr_ind_1977[22])*herb_pr$Glyph4_per_gal[22] herb_pr$Gly_pr_per_lbai_nominal <- herb_pr$Gly_pr_per_gal_nominal/4 ### all tracked prices in records refer to Glyphosate 4 lbs. a.i. per gal #### reconstruct Paraquat price history Para_08_01 <- read.csv("Para_2008_2001_prices.csv") Para_08_01 <- Para_08_01[,c("Year", "Value")] herb_pr <- merge(herb_pr, Para_08_01, by="Year", all.x=TRUE) colnames(herb_pr)[6] <- c("Para25_per_gal") ### Paraquat prices (price per gal, as "Gramoxone Extra, 2.5 lbs.ai/gal"), taken from 1991-2000 Agricultural Price Summaries scanned pdfs herb_pr[which(herb_pr$Year>=1991 & herb_pr$Year<=2000), c("Para25_per_gal")] <- c(36.10, ##1991 32.10, 32.60, 33.20, 35.00, 36.20, 37.80, 39.00, 34.80, 34.40) ##2000 ### Composite Paraquat price history using Herbicide Price Index, 1991 price as fulcrum herb_pr$Para_pr_per_gal_nominal <- herb_pr$Para25_per_gal herb_pr$Para_pr_per_gal_nominal[which(herb_pr$Year>=1970 & herb_pr$Year<=1990)] <- (herb_pr$Herb_pr_ind_1977[which(herb_pr$Year>=1970 & herb_pr$Year<=1990)]/herb_pr$Herb_pr_ind_1977[22])*herb_pr$Para25_per_gal[22] herb_pr$Para_pr_per_lbai_nominal <- herb_pr$Para_pr_per_gal_nominal/2.5 ### all tracked prices in records use Gramoxone = paraquat 2.5 lbs. a.i. per gal herb_pr[which(herb_pr$Herb_pr_ind_1977==9999),2] <- NA ### add herbicide price history to Price matrix prices <- merge(prices, herb_pr[,c("Year", "Gly_pr_per_lbai_nominal", "Para_pr_per_lbai_nominal")], by="Year", all.x=TRUE) # write.csv(prices, file="prices_final.csv") ### output nominal prices with running PPIPW composite (1970-2014) as summary .csv ### Build schedule of machinery steps ### 1) schedule of tillage years in Alternate tillage alt <- rep(c(0,0,1), length.out = 45) alt <- data.frame(cbind(seq(1970, 2014, by=1), alt)) colnames(alt) <- c("Year", "tilled") ### 2) generate list of all tillage steps taken by tillage treatment mach <- matrix(9999, nrow = dim(yld)[1], ncol=6) ### machine use activity listing order is c(Disk, MBplot, ChPlow, Cultivate, Cultimulch, HerbSpray) for (i in 1:dim(yld)[1]) { if (yld$tillage[i]==1) { ### MP mach[i,] <- c(1,1,0,1,1,0) } else if (yld$tillage[i]==2) { ### AT ifelse(alt$tilled[which(alt$Year==yld$Year2[i])]==0, mach[i,] <- c(0,0,0,0,0,1), ### in NT years mach[i,] <- c(1,1,0,1,1,0)) ### in MB years } else if (yld$tillage[i]==3) { ### Chisel till mach[i,] <- c(2,0,1,1,1,0) } else { ### NT mach[i,] <- c(0,0,0,0,0,1) } } colnames(mach) <- c("Disk", "MPlow", "ChPlow", "Cultivate", "Cultimulch", "HerbSpray") yld <- cbind(yld, mach) ### schedule of machine use steps based on fertilizer application fert <- matrix(3, nrow = dim(yld)[1], ncol=2) for (i in 1:dim(yld)[1]) { if (yld$fert[i]==1) { ### Control fert[i,] <- c(0,0) } else if (yld$fert[i]==2 | yld$fert[i]==3) { ### N-only(2) and N+NPK(3) fert[i,] <- c(1,0) } else { ### NPK(4) and NPK+NPK(5) fert[i,] <- c(1,1) } } yld <- cbind(yld, fert) colnames(yld)[22:23] <- c("Napply", "PKapply") yld[which(yld$Crop=="Soy"), c("Napply", "PKapply")] <- c(0,0) ### No fertilizer application in soybean years ### schedule of fertilizer application rates -- all rates in lbs N-P-K/ac (NOT lbs. N-P2O5-K2O) #### schedule of N application Nsched <- matrix(9999, nrow=length(unique(yld$Year2)), ncol=6) Nsched[,1] <- seq(1970, 2014, by=1) colnames(Nsched) <- c("Year", "F1", "F2", "F3", "F4", "F5") Nsched[,2] <- 0 ## Controls Nsched <- data.frame(Nsched) Nsched[which(Nsched$Year>=1970 & Nsched$Year<=1972),c("F2","F3","F4","F5")] <- 125 ## '70-72 Nsched[which(Nsched$Year>1972),c("F2","F3","F4","F5")] <- 175 ## '73-2014 Nsched[which(Nsched$Year %in% seq(1991, 2013, by=2)),2:6] <- 0 ### remove fertilizer application in soybean years #### schedle of P application, corrected to lbs P from lbs P2O5 Psched <- matrix(9999, nrow=length(unique(yld$Year2)), ncol=6) Psched[,1] <- seq(1970, 2014, by=1) colnames(Psched) <- c("Year", "F1", "F2", "F3", "F4", "F5") Psched[,2:3] <- 0 ## Controls and N-only Psched <- data.frame(Psched) Psched[which(Psched$Year>=1970 & Psched$Year<=1972),c("F3","F4","F5")] <- 50/2.29 ## '70-72 all NPK treatments Psched[which(Psched$Year>=1973 & Psched$Year<=1999),c("F3","F4","F5")] <- 80/2.29 ## '73-99 all NPK treatments Psched[which(Psched$Year>=2000),c("F3","F4","F5")] <- 50/2.29 ## '00-14 all NPK treatments Psched[which(Psched$Year %in% seq(1991, 2013, by=2)),2:6] <- 0 ### remove fertilizer application in soybean years #### schedle of K application, corrected to lbs K from lbs K2O Ksched <- matrix(9999, nrow=length(unique(yld$Year2)), ncol=6) Ksched[,1] <- seq(1970, 2014, by=1) colnames(Ksched) <- c("Year", "F1", "F2", "F3", "F4", "F5") Ksched[,2:3] <- 0 ## Controls and N-only Ksched <- data.frame(Ksched) Ksched[which(Ksched$Year>=1970 & Ksched$Year<=1972),c("F3")] <- 25/1.2 ## '70-72 N+NPK Ksched[which(Ksched$Year>=1973 & Ksched$Year<=1999),c("F3")] <- 120/1.2 ## '73-99 N+NPK Ksched[which(Ksched$Year>=1970 & Ksched$Year<=1972),c("F4","F5")] <- 125/1.2 ## '70-72 NPK/NPK+NPK Ksched[which(Ksched$Year==1973 ),c("F4","F5")] <- 240/1.2 ## '73 NPK/NPK+NPK Ksched[which(Ksched$Year>=1974 & Ksched$Year<=1999),c("F4","F5")] <- 180/1.2 ## '74-99 NPK/NPK+NPK Ksched[which(Ksched$Year>=2000),c("F3","F4","F5")] <- 150/1.2 ## '00-14, all NPK treatments Ksched[which(Ksched$Year %in% seq(1991, 2013, by=2)),2:6] <- 0 ### remove fertilizer application in soybean years ### look up fertilizer rate by year and apply to Fert Rate schedule for each treatment in each year Frate <- matrix(9999, nrow=dim(yld)[1], ncol=3) for (i in 1:dim(yld)[1]) { Frate[i,] <- c(Nsched[match(yld$Year2[i], Nsched$Year), c(paste("F",yld$fert[i], sep=""))], Psched[match(yld$Year2[i], Psched$Year), c(paste("F",yld$fert[i], sep=""))], Ksched[match(yld$Year2[i], Ksched$Year), c(paste("F",yld$fert[i], sep=""))]) } colnames(Frate) <- c("Nrate","Prate","Krate") yld <- cbind(yld, Frate) #### herbicide use schedule herb_use <- read.csv("Herb_use.csv") ## record of herbicide application in No-till herb_use[is.na(herb_use[,2]),2] <- 0 herb_use[is.na(herb_use[,3]),3] <- 0 Hrate <- matrix(9999, nrow=dim(yld)[1], ncol=2) Hrate <- as.data.frame(Hrate) for(h in 1:dim(yld)[1]){ Hrate[h,1] <- yld$HerbSpray[h]*herb_use[match(yld$Year2[h], herb_use$Year),2] Hrate[h,2] <- yld$HerbSpray[h]*herb_use[match(yld$Year2[h], herb_use$Year),3] } colnames(Hrate) <- c("Gly_lbsai", "Para_lbsai") yld <- cbind(yld, Hrate) #### Lime application history lime_app <- read.csv("Lime_appl_2014_1970.csv") ###### Machinery costs #### Machinery use $2013 cost by implement schedule (as drawn from Illinois Farmdocs) mach_op <- read.csv("Mach_op_costs.csv") ### reconstruct Machinery Prices Paid Index machind2011 <- read.csv("Mach_price_index_2011.csv") # 1990-2014 machinery total price index (2011-pegged), from USDA NASS Quick Stats machind2011 <- machind2011[, c("Year", "Value")] colnames(machind2011)[2] <- c("MACHIND11") machind1990 <- read.csv("Mach_price_index_1990.csv") # 1975-2013 machinery total price index (1990-pegged), from USDA NASS Quick Stats machind1990 <- machind1990[, c("Year", "Value")] colnames(machind1990)[2] <- c("MACHIND90") ### 1977-indexed Machinery Prices Paid Index to back calculate from 1970-1975 (from scanned pdfs of NASS annual price summaries back to 1971) machind1977 <- data.frame(c(49,## 1970 51, 54, 58, 68, 82, 91, 100, 109, 122, 136, 152, 165, 174, 181, 178, 174, 174, 181, 193, 202, 211, 219, 227)) ## 1993 machind1977$Year <- seq(1970,1993, by=1) colnames(machind1977)[1] <- c("MACHIND77") machind_comp <- merge(machind1977, machind1990, by="Year", all = TRUE) machind_comp <- merge(machind_comp, machind2011, by="Year", all=TRUE) machind_comp$MACHINDcomp <- machind_comp$MACHIND11 ### back-calculate 2011-indexed Machinery Prices Paid index to from 1989 to 1975 using the 1990-pegged mach. index (use 1990 index value as fulcrum) machind_comp$MACHINDcomp[which(machind_comp$Year<=1989 & machind_comp$Year>=1975)] <- machind_comp$MACHINDcomp[21]*(machind_comp$MACHIND90[which(machind_comp$Year<=1989 & machind_comp$Year>=1975)]/machind_comp$MACHIND90[22]) ### back-calculate 2011-indexed Machinery Prices Paid index from 1974 to 1970 using the 1977-pegged mach. index (use 1977 value as fulcrum) machind_comp$MACHINDcomp[which(machind_comp$Year<=1974 & machind_comp$Year>=1970)] <- machind_comp$MACHINDcomp[8]*(machind_comp$MACHIND77[which(machind_comp$Year<=1974 & machind_comp$Year>=1970)]/machind_comp$MACHIND77[8]) ### add composited machinery prices paid index into prices matrix prices <- merge(prices, machind_comp[,c("Year", "MACHINDcomp")], by="Year", all=TRUE) # write.csv(prices, file="prices_final.csv") #### ### compute treatment costs by year costs <- yld[,c(4:5,9,6,7,2,1,8,10:15)] ### set up costs table, calculate inflation adjustment colnames(costs)[1] <- c("Year") ### by line, multiply sum-product of activity needs by cost of each activity by machinery cost index by PPITW costs$mach_cost <- NA for(i in 1:dim(costs)[1]){ costs$mach_cost[i] <- sum(yld[i,16:23]*mach_op$Mach_op_cost_2013)*(prices$MACHINDcomp[which(prices$Year==yld$Year2[i])]/prices$MACHINDcomp[44])*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } costs$diesel_gal <- NA for(i in 1:dim(costs)[1]){ costs$diesel_gal[i] <- sum(yld[i,16:23]*mach_op$Gal.diesel) } costs$diesel_cost <- NA for(i in 1:dim(costs)[1]){ costs$diesel_cost[i] <- costs$diesel_gal[i]*prices$Diesel_pr_nominal[which(prices$Year==yld$Year2[i])]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } costs$labor_hrs <- NA for(i in 1:dim(costs)[1]){ costs$labor_hrs[i] <- sum(yld[i,16:23]*mach_op$Hrs..Labor) } costs$labor_cost <- NA for(i in 1:dim(costs)[1]){ costs$labor_cost[i] <- costs$labor_hrs[i]*prices$Wage_pr_nominal[which(prices$Year==yld$Year2[i])]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } costs$N_cost <- NA ### switched to urea after 2004 for(i in 1:dim(costs)[1]){ if(costs$Year2[i]<=2004){ costs$N_cost[i] <- yld$Nrate[i]*prices$Amm_nit_pr_nominal[which(prices$Year==yld$Year2[i])]/(0.34*2000)*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } else{ costs$N_cost[i] <- yld$Nrate[i]*prices$urea_pr_nominal[which(prices$Year==yld$Year2[i])]/(0.45*2000)*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } } costs$P_cost <- NA for(i in 1:dim(costs)[1]){ costs$P_cost[i] <- (yld$Prate[i]*prices$TSP_pr_nominal[which(prices$Year==yld$Year2[i])]*2.29)/(0.45*2000)*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } costs$K_cost <- NA for(i in 1:dim(costs)[1]){ costs$K_cost[i] <- (yld$Krate[i]*prices$KCl_pr_nominal[which(prices$Year==yld$Year2[i])]*1.2)/(0.60*2000)*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } costs$Fert_cost <- costs$N_cost+costs$P_cost+costs$K_cost costs$Herb_cost <- NA for(i in 1:dim(costs)[1]){ costs$Herb_cost[i] <- sum(yld[i, c("Gly_lbsai","Para_lbsai")]*prices[which(prices$Year==yld$Year2[i]), c("Gly_pr_per_lbai_nominal", "Para_pr_per_lbai_nominal")], na.rm=TRUE)*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==yld$Year2[i])]) } costs$Lime_cost <- NA for(p in unique(costs$Plot_Idnum)){ bulklime <- ( (lime_app$X1975[which(lime_app$Plot==p)]*prices$Lime_pr_nominal[which(prices$Year==1975)]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==1975)])) + (lime_app$X1983[which(lime_app$Plot==p)]*prices$Lime_pr_nominal[which(prices$Year==1983)]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==1983)])) + (lime_app$X1991[which(lime_app$Plot==p)]*prices$Lime_pr_nominal[which(prices$Year==1991)]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==1991)])) + (lime_app$X2001[which(lime_app$Plot==p)]*prices$Lime_pr_nominal[which(prices$Year==2001)]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==2001)])) + (lime_app$X2012[which(lime_app$Plot==p)]*prices$Lime_pr_nominal[which(prices$Year==2012)]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==2012)]))) costs$Lime_cost[which(costs$Plot_Idnum==p)] <- bulklime/length(unique(costs$Year2)) } costs$Total_overhead <- NA for(i in 1:dim(costs)[1]){ # i=1 costs$Total_overhead[i] <- sum(costs[i,c("Fert_cost", "Herb_cost", "Lime_cost", "labor_cost", "diesel_cost")]) } costs$Crop_rev <- NA for(i in 1:dim(costs)[1]){ if(costs$Crop[i]=="Corn"){ costs$Crop_rev[i] <- costs$yield_bu_ac[i]*prices$Corn_pr_nominal[which(prices$Year==costs$Year2[i])]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==costs$Year2[i])]) } else { costs$Crop_rev[i] <- costs$yield_bu_ac[i]*prices$Soy_pr_nominal[which(prices$Year==costs$Year2[i])]*(prices$PPITWcomp[45]/prices$PPITWcomp[which(prices$Year==costs$Year2[i])]) } } costs$Profit <- costs$Crop_rev-costs$Total_overhead ### calculate cumulative profits in individual treatment plots costs$cum.profit.cont <- NA ### cumulative profit in each individual treatment plot, continuous from 1970-2014 for(i in 1:dim(costs)[1]){ costs$cum.profit.cont[i] <- sum(costs$Profit[which(costs$Plot_Idnum==costs$Plot_Idnum[i] & costs$Year2<=costs$Year2[i])], na.rm=TRUE) } costs$cum.profit.rot <- NA ### cumulative profits by plot, re-set when Rotation changes for(i in 1:dim(costs)[1]){ costs$cum.profit.rot[i] <- sum(costs$Profit[which(costs$Plot_Idnum==costs$Plot_Idnum[i] & costs$Year2<=costs$Year2[i] & costs$Rotation==costs$Rotation[i])], na.rm=TRUE) } #### export results write.csv(costs, file="TxF_Econ_model_output.csv")