portfolio_6_factor/
portfolio_6_factor.rs

1
2//!
3//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
4//!
5//! File : portfolio_6_factor.rs
6//!
7//! Purpose :   Implements a portfolio optimization model using factor model.
8//!
9
10extern crate mosek;
11extern crate itertools;
12use mosek::{Task,Boundkey,Objsense,Streamtype,Soltype,Transpose,Solsta};
13use itertools::{iproduct};
14use std::convert::TryInto;
15
16
17const INF : f64 = 0.0;
18
19#[allow(non_snake_case)]
20fn portfolio(w      : f64,
21             mu     : &[f64],
22             x0     : &[f64],
23             gammas : &[f64],
24             theta  : &[f64],
25             GT     : &Matrix) -> Result<Vec<(f64,f64)>,String> {
26
27    let mut task = match Task::new() {
28        Some(e) => e,
29        None => return Err("Failed to create task".to_string()),
30        }.with_callbacks();
31    task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
32
33    let (kx,nx) = GT.size();
34    if mu.len() != nx { panic!("Mismatching data"); }
35
36    let k : i32 = kx.try_into().unwrap();
37    let n : i32 = nx.try_into().unwrap();
38    let total_budget : f64 = w + x0.iter().sum::<f64>();
39
40    //Offset of variables into the API variable.
41    let numvar = n;
42    let voff_x : i32 = 0;
43
44    // Constraint offset
45    let coff_bud : i32 = 0;
46
47    // Holding variable x of length n
48    // No other auxiliary variables are needed in this formulation
49    task.append_vars(numvar)?;
50
51    // Setting up variable x 
52    for j in 0..n {
53        task.put_var_name(voff_x+j,format!("x[{}]",j+1).as_str())?;
54    }
55    task.put_var_bound_slice_const(voff_x,voff_x+n, Boundkey::LO,0.0, INF)?;
56
57    // One linear constraint: total budget
58    task.append_cons(1)?;
59    task.put_con_name(coff_bud, "budget")?;
60
61    /* Coefficients in the first row of A */
62    for j in 0..n {
63        task.put_aij(coff_bud,voff_x + j, 1.0)?;
64    }
65    task.put_con_bound(coff_bud, Boundkey::FX, total_budget, total_budget)?;
66
67    // Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
68    // We need k+n+1 rows and we fill them in in three parts
69    task.append_afes((k+n) as i64 + 1)?;
70    // 1. The first affine expression = gamma, will be specified later
71    // 2. The next k expressions comprise G_factor_T*x, we add them row by row
72    //    transposing the matrix G_factor on the fly
73    task.put_afe_f_row_list((1..1+k as i64).collect::<Vec<i64>>().as_slice(), // f row idxs
74                            vec![n; k as usize].as_slice(), // row lengths
75                            (0..GT.len() as i64).step_by(n as usize).collect::<Vec<i64>>().as_slice(), // row ptr
76                            iproduct!(0..k,0..n).map(|(_,b)| b).collect::<Vec<i32>>().as_slice(), // varidx, 0..n repeated k times
77                            GT.data_by_row().as_slice())?;
78    // 3. The remaining n rows contain sqrt(theta) on the diagonal
79    for (i,thetai) in (0..n).zip(theta.iter()) {
80        task.put_afe_f_entry(i as i64 + 1 + k as i64, voff_x + i, thetai.sqrt())?;
81    }
82
83    // Input the affine conic constraint (gamma, GT*x) \in QCone
84    // Add the quadratic domain of dimension k+1
85    let qdom = task.append_quadratic_cone_domain(k as i64+ 1)?;
86    // Add the constraint
87    task.append_acc_seq(qdom, 0, vec![0.0; k as usize+1].as_slice())?;
88    task.put_acc_name(0, "risk")?;
89
90    // Objective: maximize expected return mu^T x
91    for (j,&muj) in (0..n).zip(mu.iter()) {
92        task.put_c_j(voff_x + j, muj)?;
93    }
94    task.put_obj_sense(Objsense::MAXIMIZE)?;
95
96    Ok(gammas.iter().filter_map(|&gamma| {
97        // Specify gamma in ACC
98        
99        task.put_afe_g(0, gamma).ok()?;
100        task.optimize().ok()?;
101        /* Display solution summary for quick inspection of results */
102        let _ = task.solution_summary(Streamtype::LOG);
103        let _ = task.write_data(format!("portfolio_6_factor-{}.ptf",gamma).as_str());
104
105        // Check if the interior point solution is an optimal point
106        if task.get_sol_sta(Soltype::ITR).ok()? != Solsta::OPTIMAL {
107            // See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
108            eprintln!("Solution not optimal!");
109            std::process::exit(1);
110        }
111
112        /* Read the results */
113        let mut xx = vec![0.0; n as usize];
114        task.get_xx_slice(Soltype::ITR, voff_x, voff_x + n, xx.as_mut_slice()).ok()?;
115        Some((gamma,xx.iter().zip(mu.iter()).map(|(&xj,&muj)| xj*muj).sum::<f64>()))
116    }).collect::<Vec<(f64,f64)>>())
117}
118
119
120
121#[allow(non_snake_case)]
122fn main() -> Result<(),String> {
123    // Since the value infinity is never used, we define
124    // 'infinity' for symbolic purposes only
125    let n : i32      = 8;
126    let w  = 1.0;
127    let mu = &[0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379];
128    let x0 = &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
129    // Factor exposure matrix, n x 2
130    let B = Matrix::new_by_row(n as usize, 2,
131                               [ 0.4256,  0.1869,
132                                 0.2413,  0.3877,
133                                 0.2235,  0.3697,
134                                 0.1503,  0.4612,
135                                 1.5325, -0.2633,
136                                 1.2741, -0.2613,
137                                 0.6939,  0.2372,
138                                 0.5425,  0.2116 ].to_vec()).unwrap();
139
140    // Factor covariance matrix, 2x2
141    let S_F = Matrix::new_by_row(2,2,
142                                 [ 0.0620, 0.0577,
143                                   0.0577, 0.0908 ].to_vec()).unwrap();
144
145    // Specific risk components
146    let theta : &[f64] = &[0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459];
147    let S_sqrt_theta = Matrix::diag_matrix(theta.iter().map(|&v| v.sqrt()).collect());
148    let P = S_F.cholesky().unwrap();
149    let BP = B.mul(&P).unwrap();
150
151    //let GT  = BP.concat_h(&S_theta.sqrt_element().unwrap()).unwrap().transpose();
152    let GT  = BP.concat_h(&S_sqrt_theta).unwrap().transpose();
153    let gammas = &[0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48];
154
155    for (gamma,expret) in portfolio(w,
156                                    mu,
157                                    x0,
158                                    gammas,
159                                    theta,
160                                    &GT)? {
161        println!("Expected return {:.4e} for gamma {:.2}", expret, gamma);
162    }
163    Ok(())
164}
165
166// Matrix with data stored in colunn format
167#[derive(Copy,Clone)]
168pub enum MatrixOrder {
169    ByRow,
170    ByCol
171}
172
173pub struct Matrix {
174    fmt  : MatrixOrder,
175    dimi : usize,
176    dimj : usize,
177    data : Vec<f64>
178}
179
180impl Matrix {
181    pub fn new(fmt : MatrixOrder, dimi : usize, dimj : usize, data : Vec<f64>) -> Option<Matrix> {
182        if dimi*dimj == data.len() {
183            Some(Matrix{fmt,dimi,dimj,data})
184        }
185        else {
186            None
187        }
188    }
189
190    pub fn new_by_col(dimi : usize, dimj : usize, data : Vec<f64>) -> Option<Matrix> {
191        Matrix::new(MatrixOrder::ByCol,dimi,dimj,data)
192    }
193    pub fn new_by_row(dimi : usize, dimj : usize, data : Vec<f64>) -> Option<Matrix> {
194        Matrix::new(MatrixOrder::ByRow,dimi,dimj,data)
195    }
196
197    pub fn size(&self) -> (usize,usize) { (self.dimi,self.dimj) }
198    pub fn len(&self) -> usize { self.data.len() }
199    pub fn diag_matrix(data : Vec<f64>) -> Matrix {
200        let mut rdata = vec![0.0; data.len()*data.len()];
201        for (r,&v) in rdata.iter_mut().step_by(data.len()+1).zip(data.iter()) { *r = v; }
202        Matrix{fmt:MatrixOrder::ByCol,dimi:data.len(),dimj:data.len(),data:rdata}
203    }
204    pub fn transpose(&self) -> Matrix {
205        Matrix {
206            fmt : match self.fmt { MatrixOrder::ByRow => MatrixOrder::ByCol, MatrixOrder::ByCol => MatrixOrder::ByRow },
207            dimi : self.dimj,
208            dimj : self.dimi,
209            data : self.data.as_slice().to_vec()
210        }
211    }
212
213    pub fn data_by_row(&self) -> Vec<f64> {
214        match self.fmt {
215            MatrixOrder::ByRow => self.data.clone(),
216            MatrixOrder::ByCol => iproduct!(0..self.dimi,0..self.dimj).map(|(i,j)| unsafe { *self.data.get_unchecked(self.dimi*j+i) }).collect()
217        }
218    }
219
220    pub fn data_by_col(&self) -> Vec<f64> {
221        match self.fmt {
222            MatrixOrder::ByCol => self.data.clone(),
223            MatrixOrder::ByRow => iproduct!(0..self.dimj,0..self.dimi).map(|(j,i)| unsafe { *self.data.get_unchecked(self.dimj*i+j) }).collect()
224        }
225    }
226
227    pub fn data_by_row_to(&self,r : &mut[f64]) {
228        match self.fmt {
229            MatrixOrder::ByRow => { r.clone_from_slice(self.data.as_slice()); },
230            MatrixOrder::ByCol => {
231                if r.len() != self.data.len() { panic!("Incorrect result array length"); }
232                for (r,(i,j)) in r.iter_mut().zip(iproduct!(0..self.dimi,0..self.dimj)) {
233                    *r = unsafe { *self.data.get_unchecked(self.dimi*j+i) }
234                }
235            }
236        }
237    }
238
239    pub fn data_by_col_to(&self,r : &mut[f64]) {
240        match self.fmt {
241            MatrixOrder::ByCol => { r.clone_from_slice(self.data.as_slice()); },
242            MatrixOrder::ByRow => {
243                for (r,(j,i)) in r.iter_mut().zip(iproduct!(0..self.dimj,0..self.dimi)) {
244                    *r = unsafe { *self.data.get_unchecked(self.dimj*i+j) };
245                }
246            }
247        }
248    }
249
250    pub fn cholesky(&self) -> Option<Matrix> {
251        if self.dimi != self.dimj { None }
252        else {
253            let mut resdata = self.data_by_col();
254            // zero data in the non-tril part
255            for ((j,i),d) in iproduct!(0..self.dimj,0..self.dimi).zip(resdata.iter_mut()) { if i < j {*d = 0.0}; }
256            if let Ok(_) = mosek::potrf(mosek::Uplo::LO,
257                                        self.dimi.try_into().unwrap(),
258                                        resdata.as_mut_slice()) {
259                Some(Matrix{fmt:MatrixOrder::ByCol,dimi:self.dimi,dimj:self.dimj,data:resdata})
260            }
261            else { None }
262        }
263    }
264
265    pub fn mul(&self,other : &Matrix) -> Option<Matrix> {
266        if self.dimj != other.dimi { None }
267        else {
268            let mut resdata = vec![0.0; self.dimi * other.dimj];
269            if let Ok(_) = mosek::gemm(match self.fmt { MatrixOrder::ByCol => Transpose::NO, MatrixOrder::ByRow => Transpose::YES },
270                                       match other.fmt { MatrixOrder::ByCol => Transpose::NO, MatrixOrder::ByRow => Transpose::YES },
271                                       self.dimi.try_into().unwrap(),
272                                       other.dimj.try_into().unwrap(),
273                                       self.dimj.try_into().unwrap(),
274                                       1.0,
275                                       self.data.as_slice(),
276                                       other.data.as_slice(),
277                                       1.0,
278                                       resdata.as_mut_slice()) {
279                Matrix::new_by_col(self.dimi,other.dimj,resdata)
280            }
281            else {
282                None
283            }
284        }
285    }
286
287    pub fn concat_h(&self, other : & Matrix) -> Option<Matrix> {
288        if self.dimi != other.dimi {
289            None
290        }
291        else {
292            let mut resdata = vec![0.0; self.data.len() + other.data.len()];
293
294            self.data_by_col_to(& mut resdata[..self.len()]);
295            other.data_by_col_to(& mut resdata[self.len()..]);
296
297            Matrix::new_by_col(self.dimi,self.dimj+other.dimj,resdata)
298        }
299    }
300}
301
302#[cfg(test)]
303mod tests {
304    #[test]
305    fn test() {
306        super::main().unwrap();
307    }
308}
portfolio_6_factor/portfolio_6_factor.rs

portfolio_6_factor/
portfolio_6_factor.rs
