$title CES production functions in MD form * CES: Constant Elasticity of Substitution * MD: Multiplicative Disturbance * Estimate coefficients in a QCP type of model set m(*) Index of data in the model, t(*) Index of observations; alias (m,mm,mm_); * Load data parameter data(t,*) Source data; $if not set in $set in in $gdxin %in%.gdx $load t < data.dim1 $load m < data.dim2 $loaddc data $gdxin * Alter information if input are zero or negative loop((t,m), abort$(data(t,m) le 0) "Negtive or zero input";); set v(m) Index of value shares; v(m) = yes$(not (ord(m) = card(m) or ord(m) = 1)); * THETA("Q"): Scale parameter, > 0 * THETA(v(m)): Value share parameters of first n-1 inputs, (0,1) * THETA("v_end"): Non-zero subsitution parameter, from -inf to 1 variable MU(t) Error terms, SSE Sum of squared errors, THETA(m) Coefficients to be estimated; equation fit_md(t) CES model in the MD model, obj Objective in the MD model; $macro agg(t) (sum(m$v(m), THETA(m)*data(t,m)**THETA("v_end")) + \ (1-sum(m$v(m), THETA(m)))*data(t,"v_end")**THETA("v_end")) obj.. SSE =e= sum(t,sqr(MU(t))); fit_md(t).. log(data(t,"Q")) =e= log(THETA("Q")) + 1/THETA("v_end")*log(agg(t)) + MU(t); model CES_MD /obj, fit_md/; * Set bounds on some unknowns THETA.LO("Q") = 0.01; THETA.UP("Q") = inf; THETA.LO(v(m)) = 0.01; THETA.UP(v(m)) = 0.99; MU.LO(t) = -1000; MU.UP(t) = 1000; * Optimize when the substitution parameter is between 0 and 1 first THETA.LO("v_end") = 0.01; THETA.UP("v_end") = 0.99; THETA.L ("v_end") = 0.5; * Save current results in "ces_md_p.gdx" CES_MD.SAVEPOINT = 1; solve CES_MD minimzing SSE using nlp; CES_MD.SAVEPOINT = 0; parameter sse_value Current value of sum of squares; sse_value = SSE.L; * Solve the problem again when the substitution parameter is negative THETA.lo("v_end") = -inf; THETA.up("v_end") = -0.01; THETA.l ("v_end") = -0.5; solve CES_MD minimzing SSE using nlp; * If the solution value is made worse, * load the previously computed solution: if (SSE.L > sse_value, execute_loadpoint 'ces_md_p.gdx';); * Report regression statistics $eval num_v card(m) - 2 set r Index of statistics report /"phi", b_1*b_%num_v%,"rho"/; parameter jacobian(t,m) Jacobian, ident(m,mm) Identity matrix, sigma2_hat Unbiased estimator for error variance, obj_md, statistics(r,*); * Get Jacobian option nlp=convertd; $echo dictMap >convertd.opt $echo jacobian >>convertd.opt CES_MD.optfile=1; solve CES_MD minimizing SSE using nlp; * Retrieve and save Jacobian from 'jacobian.gdx' $eval num_e card(t) + 1 $eval num_x card(m) + card(t) + 1 set elist Equation list /e1*e%num_e%/, xlist Variable list /x1*x%num_x%/; parameter A(*,*) Jacobian from 'jacobian.gdx', df(elist,xlist) Alias name of Jacobian; $if not defined eset sets eset(elist), xset(xlist); execute_load 'jacobian.gdx', eset=i, xset=j, df=A; set fit_md_em(elist,t) Equation mapping, theta_vm(xlist,m) Variable mapping; execute_load 'dictmap.gdx', fit_md_em, theta_vm; loop((fit_md_em(eset, t), theta_vm(xset, m)), jacobian(t,m) = df(eset, xset);); * Find the inverse of squared Jacobian through (J_t*J)*COV=sigma_hat * where J_t is the transpose of Jacobian, sigma_hat is unbiased estimator for error variance sigma2_hat = SSE.L / (card(t) - card(m)); variable DUM Dummy objective variable, INT(t,m) Intermediate results in covariance estimation, COV(m,mm) Covariance matrix; equation edummy Dummy objective function, inv1(m,mm) First step in calculating inverse of squared Jacobian, inv2(t,m) Second step in calculating inverse of squared Jacobian; edummy.. DUM =e= 0; inv1(m,mm).. sum(t, jacobian(t,m) * INT(t,mm)) =e= sigma2_hat$(ord(m) = ord(mm)) + 0$(ord(m) ne ord(mm)); inv2(t,m).. sum(mm, jacobian(t,mm) * COV(mm,m)) =e= INT(t,m); model inv /edummy,inv1,inv2/; INT.L(t,m) = 1; COV.L(m,mm) = 1; solve inv using lp minimizing DUM; * Report statistics loop((m,r)$(ord(m) = ord(r)), statistics(r,"estimator") = THETA.L(m); statistics(r,"std error") = sqrt(COV.L(m,m)); statistics(r,"T value") = THETA.L(m)/sqrt(COV.L(m,m)); * Use the BETAREG function: statistics(r,"P value") = BETAREG( (card(t) - card(m))/(card(t) - card(m) + sqr(statistics(r,"T value"))), (card(t) - card(m))/2, 0.5 ); ); * Record objective value of MD estimation obj_md = SSE.L; display statistics, obj_md;