Changeset 1406 for applications/robust

Timestamp:

11/16/11 21:02:52 (14 years ago)

Author:

sindj

Message:

Newton solution of one-step ahead control equation. JS

Files:

• applications/robust/main.cpp (modified) (5 diffs)

applications/robust/main.cpp

@@ -16 +16 @@
 #include "ddeml.h"
 #include "stdio.h"
+#include <itpp/itoptim.h>
 
 //#include "DDEClient.h"
@@ -30 +31 @@
 const double apriorno     = 0.01;
 const int max_window_size = 40;
+const int utility_order   = 19;
+const int prediction_time = 30;
 
 /*
@@ -67 +70 @@
 */
 
+double valueCRRAUtility(const double &position, const vec &samples, const int order)
+{
+        double value = 0;
+                        
+        for(int i=0;i<samples.length();i++)
+        {
+                double sample = samples.get(i);
+                value += sample/pow(position*sample+1,order+1);
+        }       
+        
+        return value;
+}
+
+double gradientCRRAUtility(const double &position, const vec &samples, const int order)
+{       
+        double value = 0;
+        
+        for(int i=0;i<samples.length();i++)
+        {
+                double sample = samples.get(i);
+                value += (-(order+1)*pow(sample,2))/pow(position*sample+1,order+2);
+        }
+
+        return value;
+}
+
+double newtonRaphson(double startingPoint, double epsilon, vec samples, int order)
+{
+        if(samples.length()>800)
+        {
+                samples.del(801,samples.size()-1);
+        }
+        
+        int count = 0;
+        
+        bool epsilon_reached = false;
+
+        while(count<1000&&!epsilon_reached)
+        {
+                double cur_value = valueCRRAUtility(startingPoint,samples,order);
+                double cur_gradient = gradientCRRAUtility(startingPoint,samples,order);
+
+                startingPoint = startingPoint - cur_value/cur_gradient;
+
+                if(cur_value<epsilon)
+                {
+                        epsilon_reached = true;
+                }
+        }
+
+        if(count==100)
+        {
+                return startingPoint; // can be different!
+        }
+        else
+        {
+                return startingPoint;
+        }
+}
+        
+        
+
+        
+
+
+
 class model
 {
@@ -307 +376 @@
         
         itpp::Laplace_RNG LapRNG = Laplace_RNG();       
-
+        
         vector<vector<string>> strings;
 
@@ -474 +543 @@
                         */              
 
-                         // PREDICTIONS
-                        pair<vec,vec> predictions = (*model_ref)->predict(3000,time,&LapRNG);
-
-                        /*
-                        cout << predictions.first << endl << endl << predictions.second << endl << "*************************************" ;
-                        pause(5);
-                        */              
-
-                        double avg_prediction = (predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()));
-
-                        (*model_ref)->predictions.ins((*model_ref)->predictions.size(),avg_prediction);
-                                        
+                        if(time>prediction_time)
+                        {                       
+                                // PREDICTIONS
+                                pair<vec,vec> predictions = (*model_ref)->predict(5000,time,&LapRNG);
+
+                                /*
+                                cout << predictions.first << endl << endl << predictions.second << endl << "*************************************" ;
+                                pause(5);
+                                */      
+
+                                double optimalInvestment = newtonRaphson(0,0.00001,predictions.second,utility_order);
+
+                                /*
+                                vec utilityValues;
+                                for(int j=0;j<1000;j++)
+                                {
+                                        utilityValues.ins(utilityValues.length(),valueCRRAUtility(-0.5+0.001*j, predictions.second, utility_order));
+                                }*/
+
+                                double avg_prediction = (predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()));
+
+                                (*model_ref)->predictions.ins((*model_ref)->predictions.size(),avg_prediction);
+                                                
+                                myfilew.open(fstring,ios::app);
+                                
+                                /*
+                                for(int j=0;j<utilityValues.length();j++)
+                                {
+                                        myfilew << utilityValues.get(j) << ",";
+                                }
+                                myfilew << endl;
+                                */
+
+                                myfilew << avg_prediction << "," << optimalInvestment << ","; 
+                                myfilew.close();
+                        }
+                        
+
+                        //preds.ins(0,data_matrix.get(0,time+1));
+                }
+
+                if(time>prediction_time)
+                {
+                        // REAL PRICE           
                         myfilew.open(fstring,ios::app);
-                        myfilew << avg_prediction << ","; 
+                        myfilew << data_matrix.get(0,time+1) << endl;
                         myfilew.close();
-                        
-
-                        //preds.ins(0,data_matrix.get(0,time+1));
-                }
-
-                
-                // REAL PRICE
-                myfilew.open(fstring,ios::app);
-                myfilew << data_matrix.get(0,time+1) << endl;
-                myfilew.close();
+                }
+