Dr. Dobb's | Metamodeling with Perl and AMPL

Metamodeling with Perl and AMPL

To create sophisticated financial models that required simulation, statistical analysis, and multistage optimization, our authors turned to AMPL and Perl.

December 01, 2004
URL:http://www.drdobbs.com/web-development/metamodeling-with-perl-and-ampl/184416171

December, 2004: Metamodeling with Perl and AMPL

Christian is president and cofounder of Elysium Digital LLC. Dessislava is an assistant professor of operations research at Babson College, MA. They can be contacted at [email protected] and [email protected], respectively.

In the 1980s, a significant conceptual breakthrough was accomplished in computational optimization with the introduction of universal optimization languages, such as the Algebraic, or Applied, Mathematical Programming Language (AMPL) and the General Algebraic Mathematical System (GAMS). AMPL and GAMS do not solve optimization problems. Instead, they provide the environment for users to describe their models in a standardized way, and call on separate solvers to solve the problem. The results are returned to human modelers via the same software.

Originally developed at AT&T Bell Labs, AMPL is a modeling language for optimization problems. With AMPL, users describe an optimization problem in the AMPL language, after which the AMPL software processes the problem and passes it to a separate optimization solver. AMPL supports numerous commercial and free solvers. For a full listing of available solvers, see http://www.ampl.com/. The Student Edition of AMPL (which limits the number of variables and constraints in the model) is included with AMPL: A Modeling Language for Mathematical Programming, Second Edition, by Robert Fourer, David M. Gay, and Brian W. Kernighan (Duxbury Press/Brooks/Cole Publishing, 2002; ISBN 0534388094). Windows, UNIX (including Linux), and Mac OS X versions of the Student Edition and compatible solvers are also available for free download from http://www.ampl.com/. The full version of AMPL can be purchased from ILOG at http://www.ilog.com/. ILOG also sells an AMPL-compatible version of CPLEX, which is one of the most powerful solvers for linear, mixed-integer, and quadratic optimization. Windows and UNIX (including Linux) binaries of AMPL and CPLEX are available. AMPL is also one of the modeling languages supported by the NEOS Server for Optimization (http://www-neos.mcs.anl.gov/neos/). The NEOS server lets users submit optimization problems via the Internet, solves them for free on a NEOS computer, and returns the results to the submitter.

A simple optimization example from the book AMPL: A Modeling Language for Mathematical Programming is as follows: A steel mill processes unfinished steel slabs into steel bands and steel coils. It can produce 200 tons of bands per hour and 140 tons of coils per hour, but cannot sell more than 6000 tons of bands per week or 4000 tons of coils per week. When sold, bands profit $25.00 per ton, and coils profit $30.00 per ton. How many tons of bands and how many tons of coils should be produced to maximize profit?

To solve this problem, you create the two input files for AMPL, like Examples 1(a) and 1(b). Of course, this is a simple problem, and AMPL can handle problems far more complex. But real-world optimization models can have infinite variety. Frequently, they require multiple refinements of the model and data, and therefore, multiple calls to optimization and data-processing software. Some optimization formulations involve solving sequences of subproblems. Other models require, in addition to optimization, the ability to analyze data through statistical techniques (to run regressions, for example), the ability to simulate scenarios, the ability to pass results from one stage of the computation as inputs for the next stages of the problem, and so on. We define such complex models as "optimization metamodels."

We ran into the problem of solving optimization metamodels while creating sophisticated financial models that required simulation, statistical analysis, and multistage optimization. As we sought to devise a good system for handling optimization metamodels, we focused on leveraging AMPL because of its modern design and its support of a wide range of solvers.

Little Languages and Optimization Systems

In a 1995 talk at Princeton University, Brian Kernighan, one of AMPL's creators, described AMPL as a prime example of a "little language." The term "little languages" was introduced for the first time by Jon Bentley (1986), who was at Bell Labs at the same time AMPL was designed. Bentley considered three different approaches to solving the problem of generating certain images:

Interactive systems (such as PhotoShop) made drawing simple images easy, but did not allow for programmatic generation of complex pictures.
Subroutine libraries worked well for generating images from big programs, but were awkward for small tasks.
Little languages provided an attractive option, since they could be powerful enough to generate complex images in an automated way, but limited enough that users could master the syntax for the language in a reasonable amount of time.

Bentley's reasoning applied perfectly to AMPL, which was designed to allow for the easy handling of both simple and complex optimization problems. We used Bentley's reasoning (as well as our own) in creating a system for optimization metamodeling.

Interactive systems were unsuitable because of the complexity of optimization metamodels. The goal of an integrated system was not to help users create simple models—this problem had already been solved by other systems. We placed a premium on supporting complexity and automation.

The little languages approach adopted by AMPL had certainly produced an effective optimization tool, and could be extended to support optimization metamodeling. This would require adding new functionality to the little language for each task that the metamodeler might need; for example, loading data, running simulations, storing and processing results in fully flexible data structures, launching external programs to process intermediate results. In fact, some such capabilities have been added recently: Today, AMPL provides not only an environment for optimization formulations but also tools for data manipulation and database access. For example, AMPL's table declaration lets you define explicit connections between sets, parameters, variables, and expressions in AMPL, and relational database tables maintained by other software, such as Microsoft Access and Excel. In addition, recently introduced command script features let you specify the number of times a model should be rerun, change values of parameters, and specify output format.

The problem with this continuous expansion of AMPL's capabilities is that a metamodel might need a great variety of functionality, to the point that the little language might not be so little anymore. If AMPL were to encompass the entire functionality of C/C++, Perl, or Java, it would be difficult to learn, to say nothing of how difficult it would be to create.

Alternatively, you could allow a little language access to outside resources by providing a mechanism for the invocation of outside code. You can imagine an AMPL program that made calls to other code written in C/C++, Perl, or Java. This option provoked the question: If one language is to be used to call code written in a second language, is it not preferable for the calling language to be the broader, general-purpose language, and the called code to be written in the special-purpose language? The two options are illustrated in Figures 1(a) and 1(b). Again, Bentley's paper shed light on the subject, in that he agrees that performance and logic suggest that flow-control elements, such as loops, are better implemented in general-purpose languages like Pascal, rather than in little languages. From a system design standpoint, we believe that the system illustrated in Figure 1(b) is the better option.

This analysis pushed us toward using a subroutine library packaged up in an API. Users can download the library, then call the functions contained in it from within his program. Inputs and outputs to an individual optimization model are variables within the main program, and can therefore be manipulated directly using the full facilities of the programming language.

An existing metamodeling system that we analyzed in this context was Matlab, a programming language that "provides core mathematics and advanced graphical tools for data analysis, visualization, and algorithm and application development" (http:// www.Mathworks.com/). One of the available toolboxes for Matlab is the Optimization Toolbox, which enables general and large-scale optimization, including linear programming, quadratic programming, nonlinear least-squares, and nonlinear equations. Accordingly, Matlab could be viewed as an example of a language that provides access to optimization functionality via a subroutine library.

However, the evolution of Matlab showed us that it is more like a little language that has been expanded to include additional functionality, instead of an efficient general-purpose programming language with subroutine library access to special functionality. This shows in the performance of the Matlab language, especially in its slow handling of loops. Since we wanted to take advantage of good performance and flexibility in our master language, we shied away from Matlab's solution.

Several candidates emerged for the programming language of our library:

Microsoft Visual Basic is easy to learn, allows for the fast creation of attractive user interfaces for programs, and has mechanisms to interact with Excel and Access. However, it has inconsistent syntax, making programming frustrating at times. Programs written in Visual Basic run slowly, and in our experience, they run unreliably, letting the same program behave differently when run repeatedly under the same conditions. This is unacceptable for many programming situations. Additionally, Visual Basic is nonportable, limiting users to Windows, and it is licensed, requiring a license for each programmer.
C/C++ was an automatic candidate, given how widely the language is used for software development. A C/C++ compiler is available on almost every platform. Programs written in C/C++ tend to run fast (unless written poorly), and there are many existing subroutine libraries available to provide programmers with additional resources. The downside of C/C++ is that programming in the language is difficult and error prone, due to the way low-level memory management and addressing are left in the hands of the programmer. Also, while there are compilers on many platforms, code usually requires some customization for each platform.
Java was a serious contender, with many strengths. The programming language is rigorously designed for large, object-oriented projects. Programs written in Java tend to be resistant to memory errors, thanks to an effective memory-management model that includes garbage collection. There are many existing Java subroutine libraries available to provide the programmer with additional resources. Finally, portability was a key goal in Java's design, and as a result, Java programs tend to run across several platforms. Java's downside relates to one of its upsides: The language's strict reliance on the object-oriented model makes learning the language difficult and writing small programs cumbersome.
Perl has very strong built-in data parsing and reformatting functionality. Programs written in Perl are portable, running on many platforms, including all major UNIX variants and Windows versions. The language allows for easy procedural programming for smaller tasks, as well as object-oriented programming for larger projects. Also, there is a wide range of add-ons already available for Perl. This was important because it told us that programmers would have many resources available to them, and that creating such an add-on ourselves was likely to be easier and well documented. On the downside, for very large projects, Perl's object-oriented model is less rigorous than that of Java. Also, graphical user interfaces are harder to create in Perl than in Visual Basic or in some C/C++ toolkits.

After some consideration, Perl emerged as the winner. Visual Basic's design weaknesses and per-programmer licensing excluded it easily. C/C++ is too cumbersome. Java was a tempting option, but the ease with which researchers can throw together Perl programs made Perl our choice.

We implemented a Perl-AMPL library and packaged it as a module that is to be distributed via the Comprehensive Perl Archive Network (CPAN; http://www.cpan.org/). CPAN, a large collection of software and documentation that add extra functionality to Perl, is complemented by the Perl Automatic Upload Server (PAUSE). This server lets creators of Perl modules upload to CPAN, where their modules become available to Perl users for free. Because PAUSE makes it possible for Perl programmers around the world to contribute to CPAN, the archive contains thousands of Perl modules on a broad range of subjects: accessing different kinds of databases, communicating via different network protocols, interoperating with different data formats, and so on.

We need to mention here that there are some commercial products that fit the main requirements of the metamodeling system-design framework we have outlined in this article. For example, ILOG's OPL Studio (http://www.ilog.com/products/oplstudio/) includes the ability to build optimization models that are then accessed from a subroutine library using Visual Basic, Java, or C/C++. The main program can generate data for the model dynamically, and the program can trigger the solving of the model and the generated data repeatedly. This indicates that OPL Studio includes the kind of functionality that we endorse for optimization metamodeling; namely, allowing the use of a general-purpose programming language that can make optimization calls to a subroutine library. However, we believe that the existence of CPAN is an important reason for selecting Perl as the general-purpose programming language of the system. By creating a framework that lets a large number of programmers contribute code that they perceive as helpful, CPAN has created the open-source framework by which a product (in this case, Perl) can outgrow the imagination of its own original creators.

The PerlAmpl Module

The PerlAmpl module was designed as an object-oriented Perl class. The AMPL object is an abstraction of an AMPL problem. The member functions let you set the problem's model, data, and options, after which member functions can be used to trigger the solver and access the results.

To use the API, include the module in your program with use Math::Ampl;, then initialize the library:

Math::Ampl::Initialize(inAmplDir, inAmplBin, inTempDir, [inPreserveFiles]);

The argument inAmplDir is the directory containing the AMPL binary, inAmplBin is the name of the binary in that directory, inTempDir is a place to store temporary files, and inPreserveFiles is an optional argument (if 1, temp files are not deleted).

You next construct an instance:

$problem = new Math::Ampl;

Then, by using $problem->Function(Arguments) syntax, you can use the following instance member functions:

Input_Mod_Append(inModText) to append inModText to the mod data, which will be passed to AMPL as the model.
Input_Mod_Clear to clear the mod data.
Input_Dat_Add_Set_Item(inSet, inDimension1, ......) adds a set item to the dat data. inSet is the name of the set. inDimension1 is the first dimension of the item, and so on.
Input_Dat_Add_Param_Item(inParam, inDimension1, ......, inValue) to add a param item to the dat data. inParam is the name of the param. inDimension1 is the first dimension of the item, and so on. inValue is the value.
Input_Dat_Clear to clear the dat data.
Input_Display_Add(inDisplayItem) to add an item to be retrieved from AMPL using Display. After solving, the values of all added items can be retrieved in Perl using Output_Display_Get_Text or Output_Display_Get_Value.
Input_Display_Clear to clear the list of Display items to retrieve.
Input_Option_Add(inOption, inValue) to add an option to be set during the Solve command, right before AMPL actually runs the solver. inOption is the option to set, inValue is the value.
Input_Option_Clear to clear the options to set during Solve.
Solve to run AMPL on the problem. It returns -1 if an error occurred, which prevented the solver from being called, and 0 if the solver was called but could not generate a solution. It returns 1 if the solver found a solution.
Solve_Best_Solver(inMinOrMax, inObjective, inSolver1, ...) to run Solve for each inSolver, keeping track of the best results. inMinOrMax must be either "min" or "max." The function uses inMinOrMax to decide if one result is better than the next. After trying all the solvers, Solve_Best_Solver reruns the best solver to fill the Display variables with the best values (unless the best solver was already the last one, in which case, it does not need to be rerun). This function has been particularly useful for solving difficult nonlinear problems where the solution found may vary from solver to solver.
Output_Display_Get_Value(inParamOrVar, inDimension1, ...) to get the retrieved value of a variable or parameter. The inParamOrVar must have been added using Input_Display_Add before calling Solve.
Output_Display_Get_Text(inParamOrVar) to get the complete text block returned by AMPL for inParamOrVar. The inParamOrVar must have been added using Input_Display_Add before calling Solve.

The PerlAmpl module addresses many of the metamodeling issues mentioned here. The Perl environment facilitates the efficient handling of data and intermediate results, allows leveraging a wide variety of modules to maximize code reuse and minimize product development time, and makes formatting end results easy. We illustrate the system's capabilities with several examples.

Example: Handling Input

While the World Wide Web contains an enormous amount of data, it is frequently difficult to download them in the appropriate format. Perl modules, such as Finance::Quote, can be used directly to download stock price data from web sites, such as Yahoo Finance. Perl modules also exist for accessing spreadsheets and databases. For example, DBD::Excel lets you access the data within an Excel file from a database interface, while Spreadsheet::TieExcel lets you tie an array within a Perl program to an Excel spreadsheet.

When a module is not available from CPAN, Perl's strong regular expression engine makes parsing unformatted (or inconveniently formatted data) an easy task. Figure 2 and Listing 1 present an example use of PerlAmpl parsing Treasury yield curve data that was then used in a computational study of credit risky bonds. The subroutine read_file reads in data on Treasuries downloaded from http://www.bloomberg.com/ (see Figure 2) and stores the yields and maturities of the 3-month, 6-month, and 1-year bills in Perl lists @gYieldList and @gPeriodList, respectively (it ignores the other information).

Example: Handling Intermediate Results and Formatting Output

We now show an example from finance—computing and plotting the mean-variance efficient frontier of a portfolio of three stocks—to illustrate how you could use the PerlAmpl module to run optimization problems multiple times, perform statistical analysis of data, plot graphs, and format output. To generate the efficient frontier, we need to solve the optimization problem in Example 2, where n is the number of stocks in the portfolio, xi are decision variables corresponding to the allocation in each of the n stocks, rtarget is the target portfolio return, and E[ri] and sij are the expected values and the covariances of asset returns i, j, i=1,...,n, j=1,...,n, respectively. sij equals the variance of return i.

Lists of expected returns, target returns, and a covariance matrix are passed to the optimization problem formulation using the PerlAmpl module. The forecasts for expected returns could be generated in Perl, for example, by running regressions or using time series techniques. A number of statistics Perl modules, such as Statistics::Regression, available free from http://www.cpan.org/, can help with the statistical analysis of data. Also, Perl modules, such as Math::Matlab, provide additional capabilities by allowing for calling outside statistical software and collecting the results.

The portfolio optimization problem is nonlinear, so the FSQP solver is called from within AMPL. The results are stored in a hash in Perl. JGraph (a free graphing software project under Linux; http://sourceforge.net/projects/jgraph/) plots the results dynamically. The advantage of this type of organization is that the whole metamodel can be run efficiently from beginning to end while keeping the same output format. Output formats that are frequently convenient, for example, are Latex table, Excel, or HTML. Listing 2 contains a subroutine for printing the optimal standard deviation for each level of target portfolio return directly in a Latex table format.

The output of this program is a file efficientFrontier.jgr and this table:

Expected Return & Standard Deviation \\ 
0.080 & 0.105 \\ 
0.085 & 0.105 \\ 
0.090 & 0.105 \\ 
0.095 & 0.107 \\ 
0.100 & 0.118 \\ 
0.105 & 0.136 \\ 
0.110 & 0.160 \\ 
0.115 & 0.186 \\ 
0.120 & 0.214 \\

The efficientFrontier.jgr file contains the text in Example 3. Using the command jgraph efficientFrontier.jgr > efficientFrontier.ps, the file can be converted to a picture (Figure 3), which can then be used, for example, in a Latex file. The conversion of the Jgraph file to a PostScript file can be automated by including it directly in the Perl program.

Example: Handling Intermediate Results And Formatting Output

We successfully used the PerlAmpl module in computational comparisons of the efficiency of multiperiod versus single-period portfolio optimization techniques. The module can be similarly used in other optimization metamodels that involve simulations at multiple stages.

Consider an example of a multiperiod optimization problem: Given a portfolio of stocks and information about their future expected values at multiple times in the future, solve two different portfolio optimization problem formulations to determine optimal portfolio allocations for each time period. Then run simulations to test which of the two allocations results in a better final period portfolio return. To compare the two optimization formulations, you would need to compare cumulative returns at the end of the time horizon. At every point in time, you need to keep track of the portfolio value, rebalance the portfolio, simulate returns for one period ahead, and recompute the value of the portfolio.

AMPL contains several built-in random generator functions, such as Beta(a,b), Cauchy(), Exponential(), Normal(mean, standard deviation), Poisson(parameter), and Uniform(lower limit, upper limit). All of the previous random distributions can also be generated in Perl by using, for example, the Math::Cephes Perl module, which contains over 150 mathematical functions. In some contexts, such as the multiperiod portfolio optimization example we mention, it is significantly more convenient to run the simulations within the main Perl program because of the need to process intermediate optimization results. Moreover, Perl offers a wider variety of random generator functions. Listing 3 is Perl code for a simplified simulation example in which the asset returns in a portfolio are drawn from a multivariate normal distribution with prespecified expected values vector and covariance matrix. Currently, there is no provision for generating correlated random variables from within AMPL. We use the Perl library random_multivariate_normal.pl (http://www.cpan.org/).

Conclusion

While the PerlAmpl module is effective and useful in its current state of development, possible enhancements include parallelizing the execution of solvers to take advantage of multiprocessors, as well as networking the PerlAmpl module to allow client computers to solve optimization metamodels by leveraging software on server computers.

TPJ

Listing 1

%kMonthHash = ("Jan" => 1, "Feb" => 2, "Mar" => 3, "Apr" => 4,  "May" => 5,
               "Jun" => 6, "Jul" => 7, "Aug" => 8, "Sep" => 9, "Oct" => 10, 
               "Nov" => 11, "Dec" => 12);
%kConvertHash = ("3month" => 90, "6month" => 180, "1year" => 360, 
               "2year" => 720, "5year" => 1800, "10year" => 3600, 
               "30year" => 10800);
sub read_file
{
    my($inFile) = @_;
    my($state, $temp, $period, $yield);
    open(INPUT, $inFile) || die("Cannot open $inFile");
    $state = 0;
    while (<INPUT>)
    {
        if ($state == 0)
        {
            #if (/U.S. Treasury yield curve/)
            if (/\*U.S. Treasuries\*/)
            {
                $state = 1;
            }
        }
        elsif ($state == 1)
        {
            #READ IN CURRENT DATE
            # Sun, 4 Jul 1999, 11:32am EDT
            if (/[A-Z][a-z][a-z],\s*([0-9]+)\s*([A-Z][a-z][a-z])\s*([0-9]+),/)
            {
              $gDay = $1;
              #assign a numerical value to the month according to %kMonthHash 
              $gMonth = $kMonthHash{$2};
              $gYear = $3;
              $state = 2;
            }
            else
            {
              die("cannot parse date");
            }
        }
        elsif ($state == 2)
        {
            #if (/Bills/)
            if (/Prc Chg/)
            {
                $state = 3;
            }
        }
        elsif (($state == 3) || ($state == 5))
        {
            if ((/Notes/) || (/\[U.S. Treasury Yield Curve\]/))
            {
                $state++;
            }
            elsif (/^(\S+)\s+[^(]+\(([^)]+)\)/) #READ IN RELEVANT DATA
            {
                $period = $kConvertHash{$1};
                push(@gPeriodList, $period);
                $yield = $2;
                push(@gYieldList, $yield);
            }
            elsif (/\s*[+-][0-9+-]+\s*$/)
            {
                # do nothing
            }
            else
            {
                die("cannot parse data: $_");
            }
        }
        elsif ($state == 4)
        {
            #if (/Bonds\s+Coupon\s+Mat\s+Date/)
            if (/Prc Chg/)
            {
                $state = 5;
            }
        }
    }
    close(INPUT);
}

Back to article

Listing 2

use Math::Ampl;
use strict;

#SPECIFICATION OF INPUT DATA:
my($NumStocks) = 3;
#vector of expected returns
my(@ExpectedReturnsList) = (0.08, 0.09, 0.12);
#vector of standard deviations for each stock
my(@StdDeviationsList) = (0.15, 0.20, 0.22);
#vector of target portfolio returns 
my(@TargetReturnsList) = (0.08, 0.085, 0.09, 0.095, 0.10, 0.105, 
                                              0.11, 0.115, 0.12);
#RESULTS STORAGE:
#hash table to store optimal portfolio standard deviation results after 
#solving all optimization problems
my(%OptimalStdDeviationsHash) = (); 
#hash table to store optimal portfolio holdings after solving the 
#portfolio optimization problem for each value of TargetReturn
my(%OptimalHoldingsHash) = (); 

#OUTPUT FILES:
my($kOutputFile) = "efficientFrontier.out"; #to store Latex table 
my($kFrontierGraphFile) = "efficientFrontier.jgr"; #to store graph

#DECLARATION OF AMPL PROBLEM INSTANCE:
my($gProblem);

#optimization model problem to be solved: minimize portfolio variance 
#subject to constraints on portfolio expected return

#optimization model file to be pased to AMPL using PerlAmpl

sub setup_ampl_mod
{
    my($mod);
    $mod =<<EODATA;
    
    param NumStocks;
    param ExpectedReturns{1..NumStocks};
    param StdDeviations{1..NumStocks};
    param TargetReturn;
    
    var holdings{1..NumStocks};

    minimize portfolio_variance:
    sum{i in 1..NumStocks}
             (StdDeviations[i]*StdDeviations[i]*holdings[i]*holdings[i]);

    subject to portfolio_target_return:
    sum{i in 1..NumStocks} (ExpectedReturns[i]*holdings[i]) >= TargetReturn;

    subject to portfolio_total:
    sum{i in 1..NumStocks}(holdings[i]) = 1;
EODATA
        $gProblem->Input_Mod_Append($mod);
}
#optimization data file to be passed to AMPL using PerlAmpl
sub setup_ampl_dat
{
    #target portfolio return passed for this instance of the problem
    my($inTargetReturn) = @_; 
    my($iStock);

    $gProblem->Input_Dat_Add_Param_Item("NumStocks", $NumStocks);
    $gProblem->Input_Dat_Add_Param_Item("TargetReturn", $inTargetReturn);
    for($iStock = 1; $iStock <= $NumStocks; $iStock++)
    {
        $gProblem->Input_Dat_Add_Param_Item("ExpectedReturns", 
                                    $iStock,$ExpectedReturnsList[$iStock-1]);
        $gProblem->Input_Dat_Add_Param_Item("StdDeviations", 
                                     $iStock,$StdDeviationsList[$iStock-1]);
    }
}
#request for AMPL to keep track of variables of interest

sub setup_ampl_display
{
    $gProblem->Input_Display_Clear;
    $gProblem->Input_Display_Add("solve_result");
    $gProblem->Input_Display_Add("portfolio_variance");
    $gProblem->Input_Display_Add("holdings");
}
#script for running the problem in AMPL and obtaining the results
sub solve_problem
{
    #target portfolio return passed for this instance of the problem
    my($inTargetReturn) = @_;
    my($solved);

    $gProblem->Input_Mod_Clear; 
    $gProblem->Input_Dat_Clear;
    setup_ampl_dat($inTargetReturn);
    setup_ampl_mod();
    setup_ampl_display();
    $solved = $gProblem->Solve;
    return $solved;
}
#PRINT GRAPH OF EFFICIENT FRONTIER DYNAMICALLY
#the output is a .jgr file which can then be converted to a .ps file

sub print_graph
{
    my($graph);
    my($targetReturn); #portfolio target return, read from list
    my($portfolioVariance); #optimal result from optimization problem
    my($portfolioStdDeviation); #to be computed from portfolio variance
    my(@CurrentHoldings); #obtained from AMPL output  
    my($iStock); #counter

    open(GRAPH, ">$kFrontierGraphFile") || 
        die ("Cannot open file \"$kFrontierGraphFile\" for graph: $!");
    $graph = "newcurve marktype none linetype solid linethickness 1";
    $graph .= " label : Efficient Frontier\n"; 
    $graph .= "\tpts ";
    foreach $targetReturn (@TargetReturnsList)
    {
        solve_problem($targetReturn);
        $portfolioVariance = 
            $gProblem-> Output_Display_Get_Value("portfolio_variance");
        $portfolioStdDeviation = sqrt($portfolioVariance);
        $graph .= " $portfolioStdDeviation $targetReturn ";
#store optimal standard deviations if necessary
        $OptimalStdDeviationsHash{"$targetReturn"} = $portfolioStdDeviation;
#if necessary, store also the optimal holdings for each value of TargetReturn
        for($iStock = 1; $iStock <= $NumStocks; $iStock++)
        {
            $CurrentHoldings[$iStock-1] =
                $gProblem->Output_Display_Get_Value("holdings",$iStock);
            $OptimalHoldingsHash{"$targetReturn,$iStock"} = 
                $CurrentHoldings[$iStock-1];
        }
    }
    print GRAPH "newgraph\n xaxis\n label : Standard Deviation\n";
    print GRAPH "yaxis \n ";
    print GRAPH "label : Expected Portfolio Return\n";
    print GRAPH "$graph\n\n";
    print GRAPH "legend defaults\n x 0.17 y 0.10\n";
    close GRAPH;
}
#PRINT A TABLE WITH RESULTS IN LATEX TABLE FORMAT
sub print_table
{
    my($portfolioStdDeviation);
    my($targetReturn);

    open(OUTPUT, ">$kOutputFile") || 
        die ("Cannot open file \"$kOutputFile\" with results: $!"); 
    printf OUTPUT "%s %s ", "Expected Return", "&";
    printf OUTPUT "%s %s \n", "Standard Deviation", "\\\\";

    foreach $targetReturn (@TargetReturnsList)
    {
        $portfolioStdDeviation = $OptimalStdDeviationsHash{$targetReturn};
        printf OUTPUT "%2.3f %s ", $targetReturn, "&";
        printf OUTPUT "%2.3f %s \n",  $portfolioStdDeviation, "\\\\";   
    }
    close OUTPUT;
}
#MAIN
#initialization
Math::Ampl::Initialize($kAmplDir, $kAmplBin, $kTempDir, 1);
$gProblem =  new Math::Ampl;

print_graph();
print_table();

Back to article

Listing 3

#Generates one path, and returns a List of Lists indexed by [time
#period, asset number]. The entries equal single-period returns for
#each stock. Single-period returns are multivariate normal random variables.

sub create_scenario
{
    #pass number of stocks in portfolio and number of time periods ahead
    my($inNumStocks, $inNumPeriods) = @_; 

    my(@SimulatedReturnsList); #single time period simulated returns
    my(@ScenarioLoL); #list of lists of asset returns for each time period
    my($iT); #time period counter 

    for ($iT = 0; $iT < $inNumPeriods; $iT++)
    {
        @SimulatedReturnsList = 
                          random_multivariate_normal($inNumStocks,
                               @ExpectedReturnsList, @CovarianceMatrixList);
        #add simulated returns for time period iT to scenario path
        @ScenarioLoL[$iT] = [ @SimulatedReturnsList ];
    }
    return (@ScenarioLoL);
}

Back to article