Sure! Thank you for responding! I have pasted my optmodel code below with my comments.
proc optmodel Printlevel=0;
/* Specifiy a variable with the number of control group rating areas */
number num_RA_control = %sysfunc(countw(&listRA_c.,%str( )));
number num_RA_all = num_RA_control+1;
/* Sets to read in data */
set<num> NOBS ; /* number of observations (predictors) */
set<num> N_RAs = 1..num_RA_control ; /* Number of control rating areas */
set<num> N_RAs_tot = 1..num_RA_all ; /* Number of control rating areas + 1 (for the treated area) */
/* Row names */
string row_names {NOBS};
/* Standardized */
/* Control group matrix */
number control_matrix {NOBS,N_RAs};
read data std_matrix_notreat_B
into NOBS=[_N_]
row_names = _NAME_
%macro read_control;
%do overRAs = 1 %to %sysfunc(countw(&listRA_c.,%str( )));
control_matrix[_N_,&overRAs.] = %scan(&listRA_c.,&overRAs.,%str( ))
%end;
;
%mend read_control;
%read_control;
/* print row_names control_matrix; */
/* Treated group vector */
number treated_vector {NOBS};
read data std_matrix_treat_B
into NOBS=[_N_]
treated_vector = /*%scan(&listRA_t.,&i.,%str( ))*/RA_FL01
;
/* print row_names treated_vector; */
/* Treated + Controls Matrix (for setting the initial guess of V) */
number matrix_all_bronze {NOBS,N_RAs_tot};
for {i in NOBS} matrix_all_bronze[i,1]=treated_vector[i];
for {i in NOBS, j in 2..num_RA_all}
matrix_all_bronze[i,j]=control_matrix[i,j-1];
/* Create initial vectors for weights (W - m rating areas X 1) and the weighting matrix (V - n predictors x n predictors) */
/* initialize all weights to be equal and sum to 1 */
number init_vec_W {N_RAs} init 1/num_RA_control;
/* std dev of the predictors */
number init_vec_V {NOBS};
for {i in NOBS} init_vec_V[i] = ((sum {j in N_RAs_tot} (matrix_all_bronze[i,j]-((sum {k in N_RAs_tot} matrix_all_bronze[i,k])/num_RA_all))**2)/(num_RA_all-1))**0.5;
for {i in NOBS} init_vec_V[i] = init_vec_V[i]/init_vec_V[9];
/* For any given V we can calculate W - This is the easy part of the problem */
number Wsol {i in N_RAs,j in 1..3} init 0; /* keep track of different solutions for W */
/* two approaches */
/* 1. Use only the MSE (assumes V=1 for all v) */
var parm_W {i in N_RAs} init init_vec_W[i] >=0 <=1;
con w_con1 : sum {i in N_RAs} parm_W[i] = 1;
min MSE = sum {i in NOBS} ((treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j]))**2;
problem v1 include parm_W MSE w_con1;
use problem v1;
solve with qp;
for {i in N_RAs,j in 1..2} Wsol[i,1]=parm_W[i];
/* 2. Use the full objective function (assumes V=init_vec_V) */
number V {i in NOBS} = init_vec_V[i];
min obj_W = sum {i in NOBS} (treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j])*(V[i])*(treated_vector[i] - sum {k in N_RAs} control_matrix[i,k]*parm_W[k]);
problem v2 include parm_W obj_W w_con1;
solve with qp;
for {i in N_RAs,j in 1..2} Wsol[i,2]=parm_W[i];
/* From the above problems you can see that W is a function of V as W(v1)~=W(v2) of the two approaches above v2 is better (obj_W<MSE) */
/*
The above methods do not get me close to getting an estimate for V.
To solve for V I need to minimize the the MSE objective function with respect to V. As an intermediate step, for any iteration of a V value we can calculate W (as done above)
Then apply that W(V) to the MSE objective function. I want to find the V that minimizes the objective MSE where W(V) is calculated by minimizing obj_W.
For example: In IML I can specify a function that I then call to minimize:
start synth(v2) global(X1_SF,X0_SF,Y1_ESI_SF,Y0_ESI_SF,F1_ESI_SF,F0_ESI_SF,X_ns_SF,Y_ESI_ns_SF);
D = 1//v2`;
V = diag(D);
** Quadratic expression **
H = F0_ESI_SF`*V*F0_ESI_SF;
** Linear expression **
f = -F1_ESI_SF`*V*F0_ESI_SF;
blc1 = j(1,ncol(F0_ESI_SF),0)||.||.;
blc2 = j(1,ncol(F0_ESI_SF),1)||.||.;
blc3 = j(1,ncol(F0_ESI_SF),1)||0||1;
blc = blc1//blc2//blc3;
init_0 = j(1,ncol(F0_ESI_SF),1/ncol(F0_ESI_SF));
optn1 = {0,0};
call NLPQUA(rc,W,H,init_0,optn1,blc,,,,f);
e = F1_ESI_SF- F0_ESI_SF*W`;
error = sum(e##2);
return(error);
AD = X_ns_SF//Y_ESI_ns_SF;
s = sqrt(var(AD`));
v_start = s/s[1];
v_start = v_start[2:ncol(v_start)];
optn = {0,0};
con1 = j(1,nrow(v_start),0);
con2 = j(1,nrow(v_start),.);
con = con1//con2;
call nlpnra(rc,result_ESI_SF,"synth",v_start,optn,con);
finish;
I am looking to do the exact same thing in OR.
*/
/* I tried writing this as a subproblem but I could not figure out how to make W a function of V */
/*
My failed attempt below (because I just cannot figure out how to make W a function of V)
NOTE: Iterating between the objective function obj_W where we hold V fixed then hold W fixed is not a way to solve this problem as V will be pushed to 0 for all
values.
*/
num epsilon = 1E-6;
num iter init 0;
num exit_crit;
num constant_V {NOBS};
var parm_V {i in NOBS} init init_vec_V[i] >=0;
min obj_Wx = sum {i in NOBS} (treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j])*(parm_V[i])*(treated_vector[i] - sum {k in N_RAs} control_matrix[i,k]*parm_W[k]);
min MSEx = sum {i in NOBS} ((treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j]))**2;
problem Wx include parm_W parm_V obj_Wx w_con1;
problem Vx include parm_W parm_V MSEx;
/* main loop */
do while (1);
/* Iterations */
iter = iter+1;
/* Keep track of V */
for {i in NOBS} constant_V[i] = parm_V[i];
/* Solve for RA weights */
use problem Wx;
fix parm_V;
solve with nlp;/* optmodel does not think this is quadratic with V fixed */
/* Solve for V (how do I make the W from problem W a function of V for this problem?) */
use problem Vx;
unfix parm_V;
fix parm_W;
solve with nlp;
unfix parm_W;
/* I would want this to calculate V - not W, but I cannot figure that out. */
/* Lets assume I got it to calculate V, then I would do this: */
/* Set an exit criteria such that when V is not changing much we exit (and calculate the final W) */
exit_crit = sum {i in NOBS} (parm_V[i]-constant_V[i])**2;
print "Iteration:" iter "Exit criteria value:" exit_crit "Objective value W:" obj_Wx "Objective value MSE:" MSEx;
if (exit_crit <= epsilon) or (iter>=50) then do;
print "Estimation complete at iteration " iter;
use problem Wx;
fix parm_V;
solve with nlp;
print "Objective value:" obj_Wx;
leave;
end;
end;
for {i in N_RAs,j in 1..3} Wsol[i,3]=parm_W[i];
num synthetic_control {NOBS};
num treat_v_synth_dif {NOBS};
for {i in NOBS} synthetic_control[i] = (sum {j in N_RAs} control_matrix[i,j]*parm_W[j]);
for {i in NOBS} treat_v_synth_dif[i] = (treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j]);
print row_names treated_vector synthetic_control treat_v_synth_dif;
/* parm_V and parm_W are our estimates - we can now apply the parm_W to all of the data to generate and outcome */
set<num> NOBS_allyrs; /* number of observations on control data */
number prems_all_years_control {NOBS_allyrs,N_RAs};
string row_names_allyrs {NOBS_allyrs};
read data prem_forOR_allyears_control_B
into NOBS_allyrs=[_N_]
row_names_allyrs = _NAME_
%macro read_control;
%do overRAs = 1 %to %sysfunc(countw(&listRA_c.,%str( )));
prems_all_years_control[_N_,&overRAs.] = %scan(&listRA_c.,&overRAs.,%str( ))
%end;
;
%mend read_control;
%read_control;
number prems_all_years_treat {NOBS_allyrs};
read data prem_forOR_allyears_treat_B
into NOBS_allyrs=[_N_]
prems_all_years_treat = RA_FL01
;
num synthetic_control_prems {NOBS_allyrs};
num treat_v_synth_dif_prems {NOBS_allyrs};
for {i in NOBS_allyrs} synthetic_control_prems[i] = (sum {j in N_RAs} prems_all_years_control[i,j]*parm_W[j]);
for {i in NOBS_allyrs} treat_v_synth_dif_prems[i] = (prems_all_years_treat[i] - sum {j in N_RAs} prems_all_years_control[i,j]*parm_W[j]);
print row_names_allyrs prems_all_years_treat synthetic_control_prems treat_v_synth_dif_prems;
quit;
