* -------------------------------------------------------------------- * Sec 1: Declarations and data input for Linear Model * -------------------------------------------------------------------- * --- SET declarations SET J Production activities / wheat, barley, rapeseed, sugarbeet /; SET I Fixed inputs / Land,Labour/; * --- Data input PARAMETER p_B(I) Fixed resource availability / Land 200 Labour 10000 / ; PARAMETER p_C(J) Cost of the different production activities / wheat 40 barley 60 rapeseed 35 sugarbeet 80 /; PARAMETER p_P(J) Prices of the different products / wheat 293 barley 503 rapeseed 319 sugarbeet 596/; PARAMETER p_R(J) Revenue of the different products ; TABLE p_A(I,J) Matrix of resource coefficients wheat barley rapeseed sugarbeet Land 1 1 1 1 Labour 25 36 27 87 ; set s Simulation steps / s1*s4/; set modes Differnt pmp alternatives /LP,Calibration, Standard, Average_cost, Paris, Elasticities, CAPRI/; set tests(modes) Differnt pmp alternatives /LP, Standard, Average_cost, Paris, Elasticities, CAPRI/; parameter p_Report(*,*,*) Sensitivity analysis; parameter p_QA_H(J) Parameter for quadratic function; parameter p_QB_H(J) Parameter for quadratic function; parameter p_lambda(J), DV(I); variable v_X(J) activities; parameter p_X_levl(modes,J) activities; parameter p_XO(J) Observed activity levels / wheat 55 barley 30 rapeseed 85 sugarbeet 30 /; * -------------------------------------------------------------------- * --- Calibration binding to reveal dual values of preferable activities * -------------------------------------------------------------------- variable v_GMcal; positive variable v_X; equations e_OBJE_H Objective function e_RESOURCE Fixed resource availability e_CALIBRATION Calibration binding; e_OBJE_H.. v_GMcal =E= sum(J,(p_P(J)-p_C(J))*v_X(J)); e_RESOURCE(I).. SUM( J, v_X(J) * p_A(I,J)) =L= p_B(I); e_CALIBRATION(J).. v_X(J) - 0.001 =L= p_XO(J); model LP /e_OBJE_H, e_RESOURCE/; solve LP using lp maximizing v_GMcal; p_X_levl("LP",J)= v_X.l(J); model reveal_dual_values /e_OBJE_H, e_RESOURCE, e_CALIBRATION/; solve reveal_dual_values using lp maximizing v_GMcal; p_lambda(J) = e_CALIBRATION.m(J); p_X_levl("Calibration",J)= v_X.l(J); * -------------------------------------------------------------------- * --- Standard - linear term of cost function equals variable costs * -------------------------------------------------------------------- variable v_GMStandard; equations e_OBJE_HH Standard approach with calibrated cost function; p_QB_H(J) = p_lambda(J)/p_XO(J); p_QA_H(J) = p_C(J); e_OBJE_HH.. v_GMStandard =E= sum(J,(p_P(J)*v_X(J)- (p_QA_H(J)+0.5*p_QB_H(J)*v_X(J))*v_X(J))); model Standard /e_OBJE_HH, e_RESOURCE/; solve Standard using nlp maximizing v_GMStandard; p_X_levl("Standard",J)= v_X.l(J); * -------------------------------------------------------------------- * --- Average cost approach * -------------------------------------------------------------------- variable v_GMAverage; parameters p_QA_HH(J), p_QB_HH(J); equations e_OBJE_HHH Howitt average cost; p_QB_HH(J) = 2*p_lambda(J)/p_XO(J); p_QA_HH(J) = p_C(J) - p_lambda(J); e_OBJE_HHH.. v_GMAverage =E= sum(J,(p_P(J)*v_X(J)- (p_QA_HH(J)+0.5*p_QB_HH(J)*v_X(J))*v_X(J))); model Average_cost /e_OBJE_HHH, e_RESOURCE/; solve Average_cost using nlp maximizing v_GMAverage; p_X_levl("Average_cost",J)= v_X.l(J); * -------------------------------------------------------------------- * --- Paris approach - no linear cost function term * -------------------------------------------------------------------- variable v_GMSParis; parameters p_QB_P(J); equations e_OBJE_P Paris cost function; p_QB_P(J)= (p_C(J)+p_lambda(J))/p_XO(J); e_OBJE_P.. v_GMSParis =E= sum(J,p_P(J)*v_X(J)- 0.5*p_QB_P(J)*v_X(J)*v_X(J)); model Paris /e_OBJE_P, e_RESOURCE/; solve Paris using nlp maximizing v_GMSParis; p_X_levl("Paris",J)= v_X.l(J); * -------------------------------------------------------------------- * --- Using elasticities - with and without predefined dual values * -------------------------------------------------------------------- parameters p_QB_E(J), p_QA_E(J); variable v_GMSelas; parameter p_elas(J) Supply elasticities / wheat 1.2 barley 1.5 rapeseed 0.7 sugarbeet 0.8 /; equations e_OBJE_E using elasticities; p_QB_E(J) = 1/p_elas(J)* p_P(J)/p_XO(J) ; p_QA_E(J) = p_C(J) + p_lambda(J)-p_QB_E(J)*p_XO(J); e_OBJE_E.. v_GMSelas =E= sum(J, p_P(J)*v_X(J) -p_QA_E(J)*v_X(J)-0.5*p_QB_E(J)*v_X(J)*v_X(J)); model Elasticities /e_OBJE_E, e_RESOURCE/; solve Elasticities using nlp maximizing v_GMSelas; p_X_levl("Elasticities",J)= v_X.l(J); * -------------------------------------------------------------------- * --- CAPRI way, with elasticities and predefined resource rents * -------------------------------------------------------------------- parameters p_QB_C(J), p_QA_C(J); parameter p_resourceRent(I) "Externally defined resource rents"/ land 150 labour 0 /; variable v_GMSCAPRI; equations e_OBJE_C CAPRI style objective using elasticities; * --- Derive slope terms from elasticity, price and observed level p_QB_C(J) = 1/p_elas(J)* p_P(J)/p_XO(J) ; * --- Treat resources as tradable items with price = desired rent p_QA_C(J) = sum(i, p_A(I,J)*p_resourceRent(I)); e_OBJE_C.. v_GMSCAPRI =E= sum(J, (p_P(J)-p_C(J))*v_X(J) -p_QA_C(J)*v_X(J)-0.5*p_QB_C(J)*v_X(J)*v_X(J)); model CAPMOD Calibraition model /e_OBJE_C/; model CAPMODQ Simulation model /e_OBJE_C, e_RESOURCE/; * --- Fix levels before calibration to obtain duals v_x.fx(J) = p_XO(J); solve CAPMOD using nlp maximizing v_GMSCAPRI; * --- Calibrate variable costs based on lambda, then release fixation of levels p_QA_C(J) = v_x.m(J); v_x.lo(J) = 0; v_x.up(J) = inf; * --- In CAPRI, we define all resources with non-zero price to be fully exhausted * for instance in calibration of milk quota rents. * You could try it for Labour in this model. * To do that, set p_resourceRent("Labour") = 2, and activate next code line. * B(I) $ p_resourceRent(I) = sum(J, XO(J)*A(I,J)); solve CAPMODQ using nlp maximizing v_GMSCAPRI; p_X_levl("CAPRI",J)= v_X.l(J); * -------------------------------------------------------------------- * --- Simulating price changes * -------------------------------------------------------------------- p_Report("LP","s0",J) = p_X_levl("LP",J); p_Report("Standard","s0",J) = p_XO(J); p_Report("Average_cost","s0",J) = p_XO(J); p_Report("Paris","s0",J) = p_XO(J); p_Report("Elasticities","s0",J) = p_XO(J); p_Report("CAPRI","s0",J) = p_XO(J); p_Report(tests,"s0","P_wheat") = p_P("wheat"); p_Report(tests,"s0","P_barley") = p_P("barley"); p_Report(tests,"s0","P_rapeseed") = p_P("rapeseed"); p_Report(tests,"s0","P_sugarbeet") = p_P("sugarbeet"); loop (s, p_P("wheat") = p_P("wheat")+40; p_P("rapeseed") = p_P("rapeseed")+60; Standard.solprint = 2; Average_cost.solprint = 2; Paris.solprint = 2; Elasticities.solprint = 2; solve LP using nlp maximizing v_GMcal; p_Report("LP",s,j) = v_X.l(J); p_Report("LP",s,"P_wheat") = p_P("wheat"); p_Report("LP",s,"P_barley") = p_P("barley"); p_Report("LP",s,"P_rapeseed") = p_P("rapeseed"); p_Report("LP",s,"P_sugarbeet") = p_P("sugarbeet"); solve Standard using nlp maximizing v_GMStandard; p_Report("Standard",s,j) = v_X.l(J); p_Report("Standard",s,"P_wheat") = p_P("wheat"); p_Report("Standard",s,"P_barley") = p_P("barley"); p_Report("Standard",s,"P_rapeseed") = p_P("rapeseed"); p_Report("Standard",s,"P_sugarbeet") = p_P("sugarbeet"); solve Average_cost using nlp maximizing v_GMAverage; p_Report("Average_cost",s,j) = v_X.l(J); p_Report("Average_cost",s,"P_wheat") = p_P("wheat"); p_Report("Average_cost",s,"P_barley") = p_P("barley"); p_Report("Average_cost",s,"P_rapeseed") = p_P("rapeseed"); p_Report("Average_cost",s,"P_sugarbeet") = p_P("sugarbeet"); solve Paris using nlp maximizing v_GMSParis; p_Report("Paris",s,j) = v_X.l(J); p_Report("Paris",s,"P_wheat") = p_P("wheat"); p_Report("Paris",s,"P_barley") = p_P("barley"); p_Report("Paris",s,"P_rapeseed") = p_P("rapeseed"); p_Report("Paris",s,"P_sugarbeet") = p_P("sugarbeet"); p_QB_E(J) = 1/p_elas(J)* p_P(J)/p_XO(J) ; p_QA_E(J) = p_C(J) + p_lambda(J)-p_QB_E(J)*p_XO(J); solve Elasticities using nlp maximizing v_GMSelas; p_Report("Elasticities",s,j) = v_X.l(J); p_Report("Elasticities",s,"P_wheat") = p_P("wheat"); p_Report("Elasticities",s,"P_barley") = p_P("barley"); p_Report("Elasticities",s,"P_rapeseed") = p_P("rapeseed"); p_Report("Elasticities",s,"P_sugarbeet") = p_P("sugarbeet"); solve CAPMODQ using nlp maximizing v_GMSCAPRI; p_Report("CAPRI",s,j) = v_X.l(J); p_Report("CAPRI",s,"P_wheat") = p_P("wheat"); p_Report("CAPRI",s,"P_barley") = p_P("barley"); p_Report("CAPRI",s,"P_rapeseed") = p_P("rapeseed"); p_Report("CAPRI",s,"P_sugarbeet") = p_P("sugarbeet"); ); execute_unload "results.gdx" p_Report; display p_Report;