How to add a constraint in PROC IML

[Bochuan]

Hi guys,

This link below is another post that resolved my question refers to the PROC IML coding, at a very beginning stage.

https://communities.sas.com/t5/Mathematical-Optimization/from-Proc-IML-to-Proc-OPTMODEL/td-p/20592

However, as following the paper (Brandt et al., 2009) provided in the link, I would like to add a weight constraint in the codes by having a little modification on the codes provided in the previous post. The aim for this step is to restrict each weight to positive given the sum of the weight to 1.

The modified part in shown here below. However, the SAS seems does not work to add the constraint on this way.

		  w=mktb+(1/countb)*charb*x`;
                  w=max(0,w)/sum(max(0,w));
		  portret=((1+w`*retb)**onemgam)/onemgam;

Could you please help me out for the problem? Thank you very much!

The full modified codes are below:

proc iml;
USE ppp.wgtfracsas;
read all var {date} into by_vars;
unique_rows=uniqueby(by_vars,1,1:nrow(by_vars)); 
read all var {logsize logbm mom} into char;
read all var {wgtfrac} into mkt;
read all var {ret} into ret;
USE ppp.equalsas;
read all var {stocks} into stocks;
onemgam=1-5;
div = J(nrow(unique_rows),1,1/nrow(unique_rows));


start opt(x) global(mkt, char, ret, stocks, onemgam, unique_rows, by_vars, div);

do i=1 to nrow(unique_rows);
          if i=nrow(unique_rows) then index=unique_rows[i]:nrow(by_vars);
          else index=unique_rows[i]:unique_rows[i+1]-1; 
          mktb=mkt[index,];
		  retb=ret[index,];
		  charb=char[index,];
		  stocksb=stocks[i,];

		  w=mktb+(1/countb)*charb*x`;
                  w=max(0,w)/sum(max(0,w));
		  portret=((1+w`*retb)**onemgam)/onemgam;


		  port=port//portret;
end;
f=port`*div;
return(f);
finish opt;

*initial values;
x0 = j(3,1,1/3);
optn={1 3};
ter={100000};

*run Sharpe Ratio;
call NLPDD(rcsr,xtheta,"opt",x0,optn,) tc=ter ; *This is the non-linear function; *rc is the vector of weights;
thetanames={'sizetheta' 'bmtheta' 'momtheta'};
create theta1 from xtheta[colname=thetanames];
append from xtheta;
quit;
run;

The original codes in the previous post that i followed is below. The selected part is the only difference between the original codes and the codes would like to add the constraints. The file "wgtfracsas" and "equalsas" could be downloaded from the previous post.

 portret=((1+(mktb+(1/stocksb)*charb*x`)`*retb)**onemgam)/onemgam;

proc iml;
USE ppp.wgtfracsas;
read all var {date} into by_vars;
unique_rows=uniqueby(by_vars,1,1:nrow(by_vars)); 
read all var {logsize logbm mom} into char;
read all var {wgtfrac} into mkt;
read all var {ret} into ret;
USE ppp.equalsas;
read all var {stocks} into stocks;
onemgam=1-5;
div = J(nrow(unique_rows),1,1/nrow(unique_rows));


start opt(x) global(mkt, char, ret, stocks, onemgam, unique_rows, by_vars, div);

do i=1 to nrow(unique_rows);
          if i=nrow(unique_rows) then index=unique_rows[i]:nrow(by_vars);
          else index=unique_rows[i]:unique_rows[i+1]-1; 
          mktb=mkt[index,];
		  retb=ret[index,];
		  charb=char[index,];
		  stocksb=stocks[i,];


		  portret=((1+(mktb+(1/stocksb)*charb*x`)`*retb)**onemgam)/onemgam;


		  port=port//portret;
end;
f=port`*div;
return(f);
finish opt;

*initial values;
x0 = j(3,1,1/3);
optn={1 3};
ter={100000};

*run Sharpe Ratio;
call NLPDD(rcsr,xtheta,"opt",x0,optn,) tc=ter ; *This is the non-linear function; *rc is the vector of weights;
thetanames={'sizetheta' 'bmtheta' 'momtheta'};
create theta1 from xtheta[colname=thetanames];
append from xtheta;
quit;
run;

[Rick_SAS · Accepted solution]

I did not look at the complete program, but I suspect the problem is the line

 w=max(0,w)/sum(max(0,w));

This statement sets w to be the scalar value max(0, w[1], w[2], ..., w[n]);

I think you are trying to do an elementwise (conditional) truncation of the w vector.

You can either use the CHOOSE function or the elementwise maximum operator.

The first option can be coded as 

w = choose(w>0, w, 0);

w = w / sum(w);

The second option can be coded as 

w = (w<>0);

w = w / sum(w);

[Bochuan]

Thank you very much for your reply Rick.

I just tried two options that you suggested. The codes can be run without any errors.

This constraint should avoid any positive w value in the results. However, SAS returns the Theta that could still generate some negative w.

I.e. If we substitute the Theta into the equation below, some of w is negative. 

          w=mktb+(1/stocksb)*charb*x`;

I am wondering should I add the constraint to scale the value for optimization by following this blog written by you as the link below? Thank you very much.

https://blogs.sas.com/content/iml/2018/08/15/optimization-nonlinear-constraints.html

[Rick_SAS]

I don't know. If you provide some sample data, then we might be able to help.

Adding a constraint to the parameter does not sound like the right approach to me. The first thing I would do is check the objective function for correctness. Please review the "Ten steps before you run an optimization."

[Bochuan]

Thanks a lot Rick. It seems I should not add a constraint to parameter here.

I attach the sample data and relative codes here. The codes are modified based on an existing post (https://communities.sas.com/t5/Mathematical-Optimization/from-Proc-IML-to-Proc-OPTMODEL/td-p/20592?n...).

The coding file “Optimization in IML for SAS_sample codes_constraint” includes scalar steps.

The file "wgtfracsas1_constraint" shows "w" could be negative even though scalar steps are applied. I think none of the "w" should be negative if we have the scalar steps.

[Rick_SAS]

Sorry, but I don't follow your logic. Your last DATA step is defining a variable that you name w, but it is not the same w that you use during the optimization. 

To see the final w that you results from the optimization, run the following code in the PROC IML block. Put it after the NLPDD call and before the QUIT statement.

start optPRINT(x) global(mkt, char, ret, stocks, onemgam, unique_rows, by_vars, div);
do i=1 to nrow(unique_rows);
          if i=nrow(unique_rows) then index=unique_rows[i]:nrow(by_vars);
          else index=unique_rows[i]:unique_rows[i+1]-1; 
          mktb=mkt[index,];
        retb=ret[index,];
        charb=char[index,];
        stocksb=stocks[i,];
          w=mktb+(1/stocksb)*charb*x`;
        w = choose(w>0, w, 0);
        w = w / sum(w);
        print "Group = " i;
        print index w;
        portret=((1+w`*retb)**onemgam)/onemgam;
        port=port//portret;
end;
f=port`*div;
return(f);
finish;

f = optPrint(xtheta);

[Bochuan]

Thank you very much Rick for you suggestion to using CHOOSE function. It successfully resolved the problem in the constrained case. It did not work well due to my fault that I made a mistake when applying suggested codes to my full dataset.

BTW, could I save the "w" to a .sas7bdat file rather than printing it out in results?

I tried to add the "CREATE" and "APPEND" as below in your codes but it could only generate w for one group.

start optPRINT(x) global(mkt, char, ret, stocks, onemgam, unique_rows, by_vars, div);
do i=1 to nrow(unique_rows);
          if i=nrow(unique_rows) then index=unique_rows[i]:nrow(by_vars);
          else index=unique_rows[i]:unique_rows[i+1]-1; 
          mktb=mkt[index,];
        retb=ret[index,];
        charb=char[index,];
        stocksb=stocks[i,];
          w=mktb+(1/stocksb)*charb*x`;
        w = choose(w>0, w, 0);
        w = w / sum(w);
        print "Group = " i;
        print index w;
		create w var {w};
		append;
        portret=((1+w`*retb)**onemgam)/onemgam;
        port=port//portret;
end;
f=port`*div;
return(f);
finish;

f = optPrint(xtheta);

[IanWakeling]

I think the problem is that you are only saving the results for the last iteration in the loop, since running 'create' will destroy any pre-existing data set.   You need to move the create statement before the loop, so something like:

  create w var {i w};
  do i = 1 to 10;
    w = ranperm(8);
    append;
  end;

[Bochuan]

Thank you very much. It works very well!