topic Re: Assistance with PROC OPTMODEL in Mathematical Optimization, Discrete-Event Simulation, and OR

Assistance with PROC OPTMODEL

Alain38 — Tue, 22 Jan 2019 10:06:51 GMT

Hello everyone,
While still practicing with proc optmodel, I tried to build a simple portfolio where each asset contributes equally to the portfolio risk.
The risk contribution (RC) of an asset depending on its weight, the problem is endogeneous. The authors of this methodology recommend to use a sequential quadratic programming algorithm and minimize the total squared differences between the RC of all pairs of assets. Everything is wrote here after in IML (optimal weights were found using Excel Solver), but my PROC OPTMODEL does not work. Perhaps the problem comes from the fact that the optimal solution W is not directly in the objective function?

proc iml;
	MVC = {0.0225 0.015 0 -.0015, 0.015  0.04  0.012 0.008,
		   0 0.012 0.09 0.021, -.0015 0.008 0.021 0.01};
	W = {0.31337624982391 0.174869881146739 0.130532835089507 0.381221033939844}; /* Results from Excel solver */
	RC = W` # ((MVC * W`) / sqrt(W * MVC * W`)); print RC;

/* Risk Contribution of the 1st Asset */
	i = 1;
	Element1 = MVC[i,1]*W[,1] + MVC[i,2]*W[,2] + MVC[i,3]*W[,3] + MVC[i,4]*W[,4];
	Var1 = W[,1]*MVC[1,1]*W[,1] + W[,1]*MVC[1,2]*W[,2] + W[,1]*MVC[1,3]*W[,3] + W[,1]*MVC[1,4]*W[,4]
		 + W[,2]*MVC[2,1]*W[,1] + W[,2]*MVC[2,2]*W[,2] + W[,2]*MVC[2,3]*W[,3] + W[,2]*MVC[2,4]*W[,4]
		 + W[,3]*MVC[3,1]*W[,1] + W[,3]*MVC[3,2]*W[,2] + W[,3]*MVC[3,3]*W[,3] + W[,3]*MVC[3,4]*W[,4]
		 + W[,4]*MVC[4,1]*W[,1] + W[,4]*MVC[4,2]*W[,2] + W[,4]*MVC[4,3]*W[,3] + W[,4]*MVC[4,4]*W[,4];
	RC1 = W[,i] * ((Element1) / sqrt(Var1)); print RC1;

/* Objective function */
	f =  (RC[1,]-RC[2,])##2 + (RC[1,]-RC[3,])##2 + (RC[1,]-RC[4,])##2
	   + (RC[2,]-RC[1,])##2 + (RC[2,]-RC[3,])##2 + (RC[2,]-RC[4,])##2
	   + (RC[3,]-RC[1,])##2 + (RC[3,]-RC[2,])##2 + (RC[3,]-RC[4,])##2
	   + (RC[4,]-RC[1,])##2 + (RC[4,]-RC[2,])##2 + (RC[4,]-RC[3,])##2;

/* Sanity check */
	RC2 = sqrt(W * MVC * W`)/ncol(MVC);
	/* RC of each asset should be equal to portfolio risk divided by the number of assets */
quit;

proc optmodel ;
/* Declare the variable */
	var W{1..4} >= 0; /* Long-only constraint */

/* Populate the model by reading in the specific data instance */
	number MVC{1..4, 1..4} = [0.0225 0.015 0     -.0015 
							  0.015  0.04  0.012 0.008 
							  0		 0.012 0.09  0.021 
							  -.0015 0.008 0.021 0.01];

	number RC{1..4}; for {k in 1..4} RC[k] = sum{i in 1..4, j in 1..4}W[k] * ((MVC[k,j]*W[j]) / sqrt(W[i]*MVC[i,j]*W[j]));

/* Minimize the objective function (total squared differences between the RC of all pairs of assets) */
	minimize f = sum{i in 1..4}(RC[i] - RC[i])^2;

/* Subject to the following constraints */
	con BUDGET: sum{i in 1..4}W[i] = 1;

/* Starting points */
	W[1] = 0.25; W[2] = 0.25; W[3] = 0.25; W[4] = 0.25;

	solve with sqp;
	print W;
quit;

Could you please help me make it work, and how could I replace the starting points (equal weights) with something more generic like:

for {i in 1..4} W[i] = 1/_N_

Thank you very much in advance for your help,

Re: Assistance with PROC OPTMODEL

RobPratt — Tue, 22 Jan 2019 21:01:16 GMT

If you want RC to depend on the variable W, you should declare it as an IMPVAR instead of a NUM. Also, your algebraic expression for RC does not match the IML matrix multiplication. Try this instead:

impvar RC{k in 1..4} = W[k] * (sum{j in 1..4}MVC[k,j]*W[j]) / sqrt(sum{i in 1..4, j in 1..4} W[i]*MVC[i,j]*W[j]);

The objective should be:

minimize f = sum{i in 1..4, j in 1..4 diff {i}}(RC[i] - RC[j])^2;

Equivalently:

minimize f = sum{i in 1..4, j in 1..4: i ne j}(RC[i] - RC[j])^2;

To initialize to equal values, do this before the solver call:

num n = 4;
for {i in 1..n} W[i] = 1/n;

Finally, the SQP solver was retired as of SAS 9.3 in July 2011. You can just do this to invoke the NLP solver:

solve;

With these changes (and with or without initializing W), the resulting solution matches your Excel output:

SAS Output

[1]	W
1	0.31337
2	0.17488
3	0.13054
4	0.38121

Re: Assistance with PROC OPTMODEL

Alain38 — Wed, 23 Jan 2019 09:09:48 GMT

@RobPrattI can't thank you enough for your huge help, this is working perfectly fine! I apologize but I would have a very last question: by adapting this code for a by-group processing, using this thread, I obtain the following error:

NOTE: Problem generation will use 4 threads.
ERROR: Out of memory during problem generation.
NOTE: Unable to create problem instance due to previous errors.
969
970      create data Work.Weights_ERC from [&byvar i] W=W_sol;
NOTE: The data set Work.WEIGHTS_ERC has 1124 observations and 3 variables.
971  quit;
ERROR: The SAS System stopped processing this step because of insufficient memory.

Given that this is a test sample with only 3 groups, I think this is more likely due to an error in my following code that the fact this optimization problem is solved using NLP, perhaps more resource intensive that in the reference thread where this is a simpler optimization problem solved with QP:

%let byvar = grp;
proc optmodel printlevel=0;
	set OBS;
	num grp {OBS};
	set <string> Assets;
	set <num,str> GROUPS_ASSETS;
	number MVC{GROUPS_ASSETS, Assets};

	read data Work.MVC into Assets=[Name];
	read data Work.MVC into OBS=[_N_] grp;
	read data Work.MVC nomiss into GROUPS_ASSETS=[k=grp i=Name] {j in Assets} <MVC[k,i,j]=col(j)>;

	set BYSET = setof {i in OBS} &byvar.[i];
	num by;

	set OBS_BY = {i in OBS: &byvar.[i] = by};
	set Assets_BY = setof {o in OBS_BY, <(grp[o]),a> in GROUPS_ASSETS} a;
	var W{Assets_BY} >= 0;
	impvar RC{l in Assets_BY} = W[l]*(sum{j in Assets_BY}MVC[by,l,j]*W[j]) / sqrt(sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j]);

	minimize f = sum{i in Assets_BY, j in Assets_BY: i ne j}(RC[i] - RC[j])^2;
	con BUDGET: sum{i in Assets_BY}W[i] = 1;

	num W_sol {GROUPS_ASSETS};
	do by = BYSET;
		put by=;
		solve;
		for {i in Assets_BY} W_sol[by,i] = W[i].sol;
	end;

	create data Work.Weights_ERC from [&byvar i] W=W_sol;
quit;

Re: Assistance with PROC OPTMODEL

RobPratt — Wed, 23 Jan 2019 14:50:36 GMT

I have a couple of suggestions to improve the efficiency.

The denominator of RC[l] does not depend on l, so you might try the following instead:

impvar Denom = sqrt(sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j]);
impvar RC{l in Assets_BY} = W[l]*(sum{j in Assets_BY}MVC[by,l,j]*W[j]) / Denom;

The objective consists of pairs of equal summands, so you can cut the number of summands in half as follows:

minimize f = 2*sum{i in Assets_BY, j in Assets_BY: i < j}(RC[i] - RC[j])^2;

If these improvements are still not enough, please supply your data.

Re: Assistance with PROC OPTMODEL

Alain38 — Thu, 24 Jan 2019 10:02:11 GMT

Thank you so much for your help. The good news is that my code was correct, the bad news is that this optimization problem is actually very resource extensive. With your more efficient code, I can at least find the optimal weights, even if this is incredibly long.

I tested with 3 groups of 1,000 assets (the maximum number of assets in my study). Please find enclosed the covariance matrix for the 1st group (in .csv due to the upload limitation). With 2 other portfolios solved with QP this is pretty quick:

minimize Variance = sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j];
con BUDGET: sum{i in Assets_BY}W[i] = 1;

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           52.32 seconds
      cpu time            1:59.94


minimize DR = 1/2 * sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j];
con VARIANCE: sum{i in Assets_BY}W[i]*MVC[by,i,i] = 1;

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           51.74 seconds
      cpu time            2:00.30

but with this portfolio solved with NLP this is extremely long:

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           30:23.27
      cpu time            1:35:21.13

In the paper here, they compare several solvers for this portfolio. In table 2, with 500 assets, we can observe a huge difference in computation time among solvers. I also noticed in their paper that the scaling constraint:

con BUDGET: sum{i in Assets_BY}W[i] = 1;

is very time consuming, so I tried to remove it, and rescale weights afterwards, but in my case this is even worse!

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           40:49.35
      cpu time            2:16:08.05

In the log I see that there were actually more iterations without the constraint: 18, 16, 20 iterations without rescaling, and 12, 12, 12 iterations with the rescale constraint.

I am still waiting for a new SAS licence that our secretary purchased a few days ago and that I will install on a more powerful computer (i5-8400 2.8GHZ and 16GB RAM), compared to currently I5-5200U 2.2GHz and 8GB RAM. So I will run more tests and the full program on the new computer.

I mark this topic as solved since the optimization code now works, but if you have more ideas to make it even more efficient I will gladly take them!

Re: Assistance with PROC OPTMODEL

RobPratt — Thu, 24 Jan 2019 00:46:47 GMT

After consulting with some of my colleagues, I have a few more suggestions that seem to help.

Replace the impvar Denom with an explicit var and con, and square both sides to avoid the sqrt:

/*   impvar Denom = sqrt(sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j]);*/
   var Denom >= 0;
   con DenomCon: Denom^2 = sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j];

Clear the denominator:

/*   con RCCon{l in Assets_BY}: RC[l] = W[l]*(sum{j in Assets_BY}MVC[by,l,j]*W[j]) / Denom;*/
   con RCCon{l in Assets_BY}: RC[l]*Denom = W[l]*(sum{j in Assets_BY}MVC[by,l,j]*W[j]);

Now the only nonlinearities are quadratic.

Re: Assistance with PROC OPTMODEL

Alain38 — Thu, 24 Jan 2019 10:01:17 GMT

Thank you very much for your suggestions and I'm sorry I just realized that in the uploaded file I forgot to add the column "Name" with _1, _2, ..., _1000 (please find it enclosed).

So I tested the following code:

%let byvar = grp;
proc optmodel printlevel=0;
	set OBS;
	num grp {OBS};
	set <string> Assets;
	set <num,str> GROUPS_ASSETS;
	number MVC{GROUPS_ASSETS, Assets};

	read data Work.MVC into Assets=[Name];
	read data Work.MVC into OBS=[_N_] grp;
	read data Work.MVC nomiss into GROUPS_ASSETS=[k=grp i=Name] {j in Assets} <MVC[k,i,j]=col(j)>;

	set BYSET = setof {i in OBS} &byvar.[i];
	num by;

	set OBS_BY = {i in OBS: &byvar.[i] = by};
	set Assets_BY = setof {o in OBS_BY, <(grp[o]),a> in GROUPS_ASSETS} a;
	var W{Assets_BY} >= 0;
	/* impvar Denom = sqrt(sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j]); */
	var Denom >= 0;
	impvar RC{l in Assets_BY} = W[l]*(sum{j in Assets_BY}MVC[by,l,j]*W[j]) / Denom;
	con DenomCon: Denom^2 = sum{i in Assets_BY, j in Assets_BY}W[i]*MVC[by,i,j]*W[j];
	con RCCon{l in Assets_BY}: RC[l]*Denom = W[l]*(sum{j in Assets_BY}MVC[by,l,j]*W[j]);

	minimize f = 2*sum{i in Assets_BY, j in Assets_BY: i < j}(RC[i] - RC[j])^2;
	con BUDGET: sum{i in Assets_BY}W[i] = 1;

	num W_sol {GROUPS_ASSETS};
	do by = BYSET;
		put by=;
		solve;
		for {i in Assets_BY} W_sol[by,i] = W[i].sol;
	end;

	create data Work.Weights_ERC from [&byvar i] W=W_sol;
quit;

And sadly this is actually worse. It's longer with more iterations. 28, 21, 26 without the constraint:

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           1:07:00.72
      cpu time            3:38:05.89

and 17, 24 and 35 with the constraint:

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           1:03:14.33
      cpu time            3:20:49.47

While comparing results, I also noticed that they were slightly different (due to the different number of iterations I presume?).

Re: Assistance with PROC OPTMODEL

RobPratt — Thu, 24 Jan 2019 13:30:00 GMT

Sorry that I failed to mention one additional change that goes along with the introduction of the RCCon constraint. Replace the impvar RC declaration with this:

var RC{Assets_BY};

Re: Assistance with PROC OPTMODEL

Alain38 — Thu, 24 Jan 2019 14:48:37 GMT

My bad, I launched the code with the modifications this morning without thinking. This is indeed much better!

NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           11:35.57
      cpu time            36:20.17

With only 7 iterations in the 3 cases.

Thank you again!

Re: Assistance with PROC OPTMODEL

RobPratt — Thu, 24 Jan 2019 18:20:05 GMT

Glad to help. By the way, which SAS/OR version are you running, and which version are you upgrading to?

Re: Assistance with PROC OPTMODEL

Alain38 — Thu, 24 Jan 2019 19:02:28 GMT

I currently have the 14.1 and I suppose I will receive the last one since this is a single-user licence. Will it make a difference in the speed of execution?

Re: Assistance with PROC OPTMODEL

RobPratt — Thu, 24 Jan 2019 21:16:29 GMT

SAS/OR 14.1 was released in July 2015. The latest production version is SAS/OR 15.1, which was released in November 2018. Performance of the NLP solver with default options on this particular instance is similar in 14.1 versus 15.1. But I would still recommend the latest version for the new features and overall performance improvements.

Re: Assistance with PROC OPTMODEL

Alain38 — Sat, 02 Feb 2019 16:17:51 GMT

Just received my new licence and this is 9.4TS1M5 with SAS/OR 14.3... I guess I will discover the 15.1 when you'll be at version 16 or higher