New Scientific Features in 2010

The goal of this page is to summarize the new or scientific features of Scilab and its environment during the year 2010.

During the year 2010, the following versions of Scilab were released :

In brief, the following is a list of the new or updated scientific features in 2010 :

Improved computational features in Scilab 5.3.0

The following informations are extracted from the Changes.

In the Optimization module, here are the changes :

Other changes are related to the statistics:

The following script computes the event x associated with the probability 0.5.

format("v",25)
p = 0.5 
q = 1-p
a = 1 
b = 2
x = cdfbet("XY",a,b,p,q)

The results are the following, with 17 significant digits:

This shows that the accuracy of the cdfbet function is now close to the maximum.

16th February 2010 : Make Matrix

A collection of test matrices.

The goal of this toolbox is to provide a collection of test matrices. These matrices might be used to test algorithms for the computation of the solution of linear systems of equations, eigenvalues, matrix norms and other dense linear algebra problems.

The current toolbox is able to generate the following matrices: dingdong, Hadamard, inverse Hilbert, Hilbert, magic, Toeplitz, Vandermonde, Pascal, border, Cauchy, circul, diagonali, Frank, Hankel, identity, Moler, Rosser, urandom, Wilkinson -, Wilkinson +.

This module provides the following functions.

This module is available in ATOMS :

To install it, type :

atomsInstall('makematrix')

In the following example, we create a Cauchy matrix of size 5.

-->A = gallery("cauchy" , 1:5 )
 A  =
 
    0.5          0.3333333    0.25         0.2          0.1666667  
    0.3333333    0.25         0.2          0.1666667    0.1428571  
    0.25         0.2          0.1666667    0.1428571    0.125      
    0.2          0.1666667    0.1428571    0.125        0.1111111  
    0.1666667    0.1428571    0.125        0.1111111    0.1        

In the following example, we compare the estimated condition number of Hilbert's matrix for increasing values of n with its theoretical value.

nv = 1:30;
for n = nv
  c(n) = log(cond(makematrix_hilbert(n)));
  e(n) = log(exp(3.5*n));
end
scf();
plot ( nv , c , "bo" );
plot ( nv , e , "r-" );
legend(["cond" "exact"]);
xtitle("Condition number of the Hilbert matrix","N","Condition");

This produces the following figure.

makematrix_export.svg

16th February 2010 : Quapro - Linear and Linear Quadratic Programming

This toolbox defines linear quadratic programming solvers. The matrices defining the cost and constraints must be full, but the quadratic term matrix is not required to be full rank.

This toolbox already existed in previous versions of Scilab. The novelty is that it is now available as a ATOMS module.

Features:

http://atoms.scilab.org/toolboxes/quapro

To install it, type :

atomsInstall('quapro')

The linpro function can solve linear programs in general form:

Minimize c'*x
A*x   <= b
Aeq*x  = beq
lb <= x <= ub

The following example is extracted from "Operations Research: applications and algorithms", Wayne L. Winstons, Section 5.2, "The Computer and Sensitivity Analysis", in the "Degeneracy and Sensitivity Analysis" subsection. We consider the problem:

Min -6*x1 - 4*x2 - 3*x3 - 2*x4  
such as:
2*x1 + 3*x2 +   x3 + 2*  x4 <= 400
  x1 +   x2 + 2*x3 +     x4 <= 150
2*x1 +   x2 +   x3 + 0.5*x4 <= 200
3*x1 +   x2 +            x4 <= 250;
x >= 0;

The following script allows to solve the problem.

c = [-6 -4 -3 -2]';
A = [
     2 3 1 2
     1 1 2 1
     2 1 1 0.5
     3 1 0 1
     ];
b = [400 150 200 250]';
ci=[0 0 0 0]';
cs=[%inf %inf %inf %inf]';
[xopt,lagr,fopt]=linpro(c,A,b,ci,cs)

This produces :

xopt = [50,100,2.842D-14,0]

14th April 2010 : Cobyla Optimization Toolbox

COBYLA is a derivative free non linear constrained optimization method.

The author of this Scilab module is Yann Collette.

This software is a Scilab interface of COBYLA2, a contrained optimization by linear approximation package developed by Michael J. D. Powell in Fortran. The original source code can be found at:

http://plato.la.asu.edu/topics/problems/nlores.html

The toolbox is managed under ATOMS :

http://atoms.scilab.org/toolboxes/scicobyla

The sources are available at:

http://forge.scilab.org/index.php/p/scicobyla

The calling sequence of the cobyla function is:

[x_opt, eval_func] = cobyla(x0, func, nb_constr, rhobeg, rhoend, message, eval_func_max)

The following is an example from the help page.

nb_constr_test_1 = 0;
xopt_test_1 = [-1 0]';

function [f, con, info] = test_1(x)
  d__1 = x(1) + 1.;
  d__2 = x(2);
  f = d__1 * d__1 * 10. + d__2 * d__2;
  con = 0;
  info = 0;
endfunction

rhobeg = 1;
rhoend = 1e-3;
message_in = 0;
eval_func_max = 200;

x0 = ones(xopt_test_1);
[x_opt, status, eval_func] = cobyla(x0, test_1, nb_constr_test_1, rhobeg, rhoend, message_in, eval_func_max);

The previous script produces the following output.

-->[x_opt, status, eval_func] = cobyla(x0, test_1, nb_constr_test_1, rhobeg, rhoend, message_in, eval_func_max)
size_x = 2 res = 0
x_opt[0] = -1.000314
x_opt[1] = -0.000438
 eval_func  =
    59.  
 status  =
    0.  
 x_opt  =
  - 1.0003138  
  - 0.0004382  

21st April 2010 : Floating Point toolbox in ATOMS

The goal of this toolbox is to provide a collection of algorithms for floating point number management. This is more a learning tool than an operationnal component, although it might complement some features which are not provided by Scilab. The flps_systemgui function allows to see the distribution of floating point numbers graphically. By varying the rounding mode, the precision and the logscale, we can actually see the distribution of a toy floating point number. Adding or suppressing the denormals can be instructive too.

The functions allow to compute automatically the properties of the current Scilab system with respect to the doubles. It is similar in spirit to the number_properties functions, except that our functions are based on macros and that the returned values are made consistent with the references cited in the bibliography. Moreover, the flps_radix function returns the current radix and allows to check that the current rounding mode is round-to-nearest (which is IEEE's default).

The functions allow to create virtual floating point systems, which allows to see their discrete nature in simplified examples. The rounding mode of such a virtual floating point system can be configured to one of the four mode from the IEEE standard. This feature allows to see the effect of the rounding mode on the distribution of floating point numbers. This is not easy to see with a straightforward Scilab, since the rounding mode is round-to-nearest, most of the time.

The following is a list of the current functions :

This toolbox is available in ATOMS :

http://atoms.scilab.org/toolboxes/floatingpoint

and is managed in the Scilab Forge :

http://forge.scilab.org/index.php/p/floatingpoint/

In order to install it, type :

atomsInstall('floatingpoint')

The flps_radix function allows to get the radix and the rounding mode of the current Scilab system. [ radix , rounding ] = flps_radix ()

On typical systems, we get the following session, meaning that we have a base-2 machine with round-to-nearest rounding mode. In other rounding modes, we would get rounding=%f, but this has never been observed by us in practice.

-->[ radix , rounding ] = flps_radix ()
 rounding  =
  T  
 radix  =
    2.  

The flps_IEEEsingle function allows to create a virtual single precision floating point system based on the IEEE standard single precision. The flps_numberformat function returns the floating point number corresponding to the given double and the given floating point system. The following script allows to produce the single precision floating point number associated with 1/3.

flps = flps_IEEEsingle ( );
flpn = flps_numberformat ( flps , 1/3 )

The previous script produces the following output. As we can see, IEEE single precision floating point numbers are associated with radix 2, precision p=24 and exponent range from -126 to 127.

-->flps = flps_IEEEsingle ( )
 flps  =
 
Floating Point System:
======================
radix= 2
p=     24
emin=  -126
vmin=  1.175D-38
emax=  127
vmax=  1.701D+38
eps=   0.0000001
r=     1
gu=    T
alpha= 1.401D-45
 
-->flpn = flps_numberformat ( flps , 1/3 )
 flpn  =
 
Floating Point Number:
======================
s= 0
M= 11184811
m= 1.3333334
e= -2
d= [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1]
flps= floating point system
======================
x= (-1)^0 * 1.3333334 * 2^-2
x= 11184811 * 2^(-2-24+1)
x= (-1)^0 * (1.[0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1])_2 * 2^-2

The following example displays a toy system, including denormals and negative numbers.

radix = 2;
p = 3;
e = 3;
flps = flps_systemnew ( "format" , radix , p , e );
h = flps_systemgui ( flps );

This produces the following figure.

toysystem_all.svg

1st of June 2010 : Optkelley

These Scilab files are implementations of the algorithms from the book 'Iterative Methods for Optimization', published by SIAM, by C. T. Kelley. The book, which describes the algorithms, is available from SIAM (service@siam.org).

The optimization codes have the calling convention [f,g] = objective(x) returns both the objective function value f and the gradient vector g. I expect g to be a column vector. The Nelder-Mead, Hooke-Jeeves, Implicit Filtering, and MDS codes do not ask for a gradient.

This module was already available in previous versions of Scilab. An important update of the module has been performed. The help pages have been entirely updated. Several bugs have been fixed and unit tests were created, along with more demonstrations scripts.

This toolbox is available in ATOMS :

http://atoms.scilab.org/toolboxes/optkelley

To install it, type :

atomsInstall('optkelley')

This toolbox provide the following algorithms:

The following is a sample session.

-->function y = quad ( x )
-->x1 = [1 1]'
-->x2 = [2 2]'
-->y = max ( norm(x-x1) , norm(x-x2) )
-->endfunction
 
-->x0 = [-1.2 1]';
 
-->v1 = x0 + [0.1 0]';
 
-->v2 = x0 + [0 0.1]';
 
-->v0 = [x0 v1 v2];
 
-->[x,lhist,histout] = optkelley_mds(v0,quad,1.e-4,100,100)
 histout  =
 
    3.     3.2572995    0.0953114    0.1414214  
    7.     3.0675723    0.1897272    0.2        
    11.    2.6925824    0.3749899    0.5656854  
    15.    1.9723083    0.7202741    0.8        
    19.    1.0049876    0.9673207    2.2627417  
    22.    0.9219544    0.4234080    1.1313708  
    26.    0.7810250    0.3006404    0.5656854  
    30.    0.7810250    0.1409295    0.2828427  
    34.    0.7810250    0.0675032    0.1414214  
    37.    0.7211103    0.0851155    0.1        
    41.    0.7211103    0.0421066    0.05       
    45.    0.7211103    0.0209308    0.025      
    49.    0.7211103    0.0104335    0.0125     
    53.    0.7211103    0.0052086    0.00625    
    57.    0.7211103    0.0026022    0.003125   
    61.    0.7211103    0.0013006    0.0015625  
    65.    0.7211103    0.0006502    0.0007812  
    69.    0.7211103    0.0003251    0.0003906  
    73.    0.7211103    0.0001625    0.0001953  
    77.    0.7211103    0.0000813    0.0000977  
 lhist  =
 
    20.  
 x  =
 
    1.4    1.4          1.3999023  
    1.6    1.5999023    1.6        

1st June 2010 : Unconstrained Optimization Problem Toolbox

The goal of this toolbox is to provide unconstrained optimization problemsin order to test optimization algorithms.

The More, Garbow and Hillstrom collection of test functions is widely used in testing unconstrained optimization software. The code for these problemsis available in Fortran from the netlib software archives.

Features

This toolbox is available in ATOMS :

http://atoms.scilab.org/toolboxes/uncprb

and is manage under Scilab's Forge :

http://forge.scilab.org/index.php/p/uncprb/

To install it, type :

atomsInstall('uncprb')

In the following session, we compute the function and Jacobian matrix of Rosenbrock's test case.

-->nprob = 1
 nprob  =
 
    1.  
 
-->[n,m,x0]=uncprb_getinitf(nprob)
 x0  =
 
  - 1.2  
    1.   
 m  =
 
    2.  
 n  =
 
    2.  
 
-->option = 3
 option  =
 
    3.  
 
-->[fvec,J]=uncprb_getfunc(n,m,x0,nprob,option)
 J  =
 
    24.    10.  
  - 1.     0.   
 fvec  =
 
  - 4.4  
    2.2  

17th June 2010 : Low Discrepancy Sequences

The goal of this toolbox is to provide a collection of low discrepancy sequences. These random numbers are designed to be used in a Monte-Carlosimulation. For example, low discrepancy sequences provide a higherconvergence rate to the Monte-Carlo method when used in numericalintegration. The toolbox takes into account the dimension of the problem, i.e.generate vectors with arbitrary size.

The current prototype has the following features :

Overview of sequences

This module currently provides the following functions:

Provides the following functions to extend the maximum dimension of the Halton and Faure sequences

Provides the following functions to suggest expert settings for the sequences :

This component currently provides the following sequences:

This toolbox is available in ATOMS at :

http://atoms.scilab.org/toolboxes/lowdisc

and is managed in the Scilab Forge :

http://forge.scilab.org/index.php/p/lowdisc/

To install it, type :

atomsInstall('lowdisc')

The following example plots the 2D Faure sequence.

lds = lowdisc_new("fauref");
lds = lowdisc_configure(lds,"-dimension",2);
lds = lowdisc_startup (lds);
[lds,computed] = lowdisc_next (lds,100);
lds = lowdisc_destroy(lds);
plot(computed(:,1),computed(:,2),"bo");
xtitle("Faure sequence","X1","X2");

This produces the following figure.

lowdisc_faure.svg

24th of June 2010: Scilab Image and Video Processing toolbox v0.5.3

SIVP intends to do image processing and video processing tasks. SIVP is meant to be a useful, efficient, and free image and video processing toolbox for Scilab.

The authors of this module are Shiqi Yu, Jia Wu, Shulin Shang and Vincent Etienne.

This module provides the following functions.

To install this module, type:

atomsInstall('SIVP')

The following script is the example of the detectlefteyes function.

    SIVP_PATH = getSIVPpath();
    im = imread(SIVP_PATH + 'images/lena.png');
    face = detectfaces(im);
    rect = face*diag([1,1,1,0.7]);
    subface = imcrop(im, rect);
    
    leyes = detectlefteyes(subface);
    [m,n] = size(leyes);
    for i=1:m,  
        im = rectangle(im, leyes(i,:)+rect*diag([1,1,0,0]), [0,255,0]);
    end;
    imshow(im);    

This produces the following figure.

SIVP

July 2010: Mmodd for Partial Differential Equations

The MmodD project contains tools for the study and solution of partial differential equations (PDEs) in 2d and 3d. A set of command-line functions and a graphical user interface let you preprocess, solve, and postprocess generic PDEs for a broad range of engineering and science applications.

http://forgesitn.univ-lyon1.fr/projects/mmodd

The project is leaded by Thierry Clopeau, with the help of the following developpers :

The module provides the help pages of the following functions:

The following functions are available:

15th July 2010 - New module : A Toolbox for Unconstrained Global Optimization of Polynomial functions

"Many problems in science and engineering can be reduced to the problem of finding optimum bounds for the range of a multivariable polynomial on a specified domain. Local optimization is an important tool for solving polynomial problems, but there is no guarantee of global optimality. For polynomial optimization problems, an alternate approach is based on the Bernstein form of the polynomial. If a polynomial is written in the Bernstein basis over a box, then the range of the polynomial is bounded by the values of the minimum and maximum Bernstein coefficients. Global optimization based on the Bernstein form does not require the iterative evaluation of the objective function. Moreover, the coefficients of the Bernstein form are needed to be computed only once, i.e., only on the initial domain box. The Bernstein coefficients for the subdivided domain boxes can then be obtained from the initial box itself. Capturing these beautiful properties of the Bernstein polynomials, global optimum for the polynomial on the given domain can be obtained. The toolbox is developed based on the above ideas."

The authors of this module are Dhiraj B. Magare, Bhagyesh V. Patil and P. S. V. Nataraj.

This toolbox is available in ATOMS at :

http://atoms.scilab.org/toolboxes/Global_Optim_toolbox/

To install it, type :

atomsInstall('Global_Optim_toolbox')

This module provides the following functions :

30th of August 2010: Scilab Wavelet Toolbox v0.1.11

This toolbox is aimed to mimic matlab wavelet toolbox. Most of the functions are similiar to their counterparts in Matlab equivalents.

The author of this module is Holger Nahrstaedt.

To install it, type:

atomsInstall('swt')

This module provides the following functions.

7th of September 2010: ANN Toolbox

This is a toolbox for artificial neural networks.

The author of this module is Ryurick M. Hristev.

Features:

The module provides the following functions.

This toolbox has long been available in Scilab. It was provided in the Former Toolbox Center, since 2005. One origin of this toolbox is the book:

"The ANN Book", Edition 1, Ryurick M. Hristev, 1998

One other related document is the thesis submitted by Ryurick M. Hristev for the degree of Master of Science in Mathematics in the University of Canterbury:

"Matrix Techniques in Articifial Neural Networks", Ryurick M. Hristev, University of Canterbury, 2000

The toolbox provides 9 demonstrations:

The demonstration "encoder 4-3-4 on ANN without biases" produces the following output.

-->// Loose 4-3-4 encoder on a backpropagation network without biases
-->// (Note that the tight 4-2-4 encoder will not work without biases)
-->// ensure the same random starting point
-->rand('seed',0);
-->// network def.
-->//  - neurons per layer, including input
-->N  = [4,3,4];
-->// inputs
-->x = [1,0,0,0;
-->     0,1,0,0;
-->     0,0,1,0;
-->     0,0,0,1]';
-->// targets, at training stage is acts as identity network
-->t = x;
-->// learning parameter
-->lp = [8,0];
-->// init randomize weights between
-->r = [-1,1];
-->W = ann_FF_init_nb(N,r);
-->// 500 epochs are enough to ilustrate
-->T = 500;
-->W = ann_FF_Std_online_nb(x,t,N,W,lp,T);
-->// full run
-->ann_FF_run_nb(x,N,W)
 ans  =
 
    0.9797963    0.0000130    0.0149148    0.0206898  
    0.0000013    0.9782197    0.0171215    0.0172901  
    0.0156931    0.0209699    0.9786566    0.0000488  
    0.0178893    0.0183367    0.0000022    0.9758474  
-->// encoder
-->encoder = ann_FF_run_nb(x,N,W,[2,2])
 encoder  =
 
    0.9830039    0.0220232    0.9736080    0.0137882  
    0.9259607    0.2809842    0.0144560    0.9837462  
    0.0157942    0.988017     0.8082282    0.2658832  
-->// decoder
-->decoder = ann_FF_run_nb(encoder,N,W,[3,3])
 decoder  =
 
    0.9797963    0.0000130    0.0149148    0.0206898  
    0.0000013    0.9782197    0.0171215    0.0172901  
    0.0156931    0.0209699    0.9786566    0.0000488  
    0.0178893    0.0183367    0.0000022    0.9758474  

This toolbox is provided under GPL licence.

This toolbox is available on ATOMS:

http://atoms.scilab.org/toolboxes/ANN_Toolbox

To install it:

atomsInstall('ANN_Toolbox')

On the negative side, there is no example in the help pages, which are sometimes poorly formatted.

14th of September 2010: Factorization of Structured Matrices Toolbox

The Factorization of Structured Matrices Toolbox is a Scilab 5 toolbox for the fast factorization of matrices with displacement structure like, e.g., (block) Toeplitz matrices.

The author of this module is Sander Wahls.

Implemented algorithms:

regular hermitian matrices with displacement structure of Stein type

Features:

Current Limitations:

This module provides the following functions.

Here is a sample session.

-->// create indefinite hermitian Toeplitz matrix
 
-->T = toeplitz(1:5,1:5)+toeplitz(%i*(0:4),-%i*(0:4));
 
-->// compute displacement with respect to the shift matrix
 
-->F = [spzeros(1,5);speye(4,4) spzeros(4,1)];
 
-->D = T-F*T*F';
 
-->// compute generator
 
-->G = [T(:,1) [0;T(2:$,1)]];
 
-->J = [1 0;0 -1];
 
-->// test generator
 
-->disp(D-G*J*G');
 
    0    0    0    0    0  
    0    0    0    0    0  
    0    0    0    0    0  
    0    0    0    0    0  
    0    0    0    0    0  
 
-->// compute LDL factorization
 
-->[L,d] = ldl_gs(F,G,1,1)
 d  =
 
    1.  
  - 1.  
  - 1.  
  - 1.  
  - 1.  
 L  =
 
  - 1.          0                         0                         0                         0                       
  - 2. - i      1.7888544 + 0.8944272i    0                         0                         0                       
  - 3. - 2.i    2.6832816 + 1.3416408i  - 1.5491933 - 0.7745967i    0                         0                       
  - 4. - 3.i    3.5777088 + 1.7888544i  - 2.0655911 - 1.0327956i    1.4605935 + 0.7302967i    0                       
  - 5. - 4.i    4.472136 + 2.236068i    - 2.5819889 - 1.2909944i    1.8257419 + 0.9128709i  - 1.4142136 - 0.7071068i  
-->// test factorization
 
-->disp(T-L*diag(d)*L');
 
    0    0                         0                         0                         0                       
    0  - 4.441D-16                 2.220D-16 + 1.332D-15i  - 8.882D-16 - 2.220D-15i  - 7.105D-15 - 1.332D-15i  
    0    2.220D-16 - 1.332D-15i  - 3.553D-15 + 2.220D-16i  - 1.776D-15 - 4.219D-15i  - 1.954D-14 - 9.770D-15i  
    0  - 8.882D-16 + 2.220D-15i  - 1.776D-15 + 4.441D-15i  - 1.332D-15 - 4.441D-16i  - 2.087D-14 - 7.994D-15i  
    0  - 7.105D-15 + 1.332D-15i  - 1.954D-14 + 8.882D-15i  - 2.087D-14 + 7.772D-15i  - 5.596D-14 + 1.554D-15i  

7th of November 2010: Linear System Inversion Toolbox v1.0.2

The Linear System Inversion Toolbox is a Scilab 5 toolbox for stable inversion of linear time-invariant systems. The inverses are optimal in the sense that some norm criterion is minimized. Decision delays are also implemented for some cases.

The author of this module is Sander Wahls.

Changes in Version 1.0.2:

* Shorter startup message * Tests are now integrated with Scilabs testing system * Improved documentation

To install it, type:

atomsInstall('lsitbx')

This module provides the following functions:

The following is a sample session.

-->// Example 1: Approximate delay-one left inverse of H(z)=[1/z;1-1/z]
-->z = poly(0,"z");
-->H = [1/z ; 1-1/z];           // create transfer function of the channel
-->SYSH = tf2ss(H);             // convert into state-space model
-->SYSH.dt = "d";                       // mark as discrete-time
-->SYSL = tf2ss(1/z);           // create state-space model of the 
-->SYSL.dt = "d";                       // mark as discrete-time
-->Q = 1;                               // covariance of the data
-->R = 1e-8*eye(2,2);           // covariance of the noise is very small
-->                             // => filter is approximately a left inverse
-->SYSK = wkf(SYSL,SYSH,Q,R);   // compute wiener-kalman filter
-->disp(clean(ss2tf(SYSK*SYSH)));       // should be approximately L(z)=1/z
 
                                       2  
  - 2.361D-09 + 1.0000000z - 6.180D-09z   
    -----------------------------------   
                      2                   
                     z                    

3rd of November 2010: Identification Toolbox v1.0

Identification Toolbox is used to construct mathematical models of LTI systems from measured input-output data sequences. The tool allows preprocessing of signals, identification of LTI systems and validating of constructed models. The tool is aimed particularly at black-box identification.

The author if this module is Martin Novotny.

This module provides the following functions.

To install it, type:

atomsInstall('identification')

30th of October 2009: HYDROGR

HYDROGR is a set of models and function to perform data analysis and modelling of hydrological time series. In particular HYDROGR contains the rainfall-runoff model GR4J developped in CEMAGREF (see http://www.cemagref.fr/webgr/index.htm). The functions are written in Scilab languages and C.

The author of this module is Julien LERAT from CEMAGREF.

This module provides the following functions.

The following is an example for the c_GR4J function.

      // P is a vector containing daily rainfall values
      P = max(0,exp(convol(1/3*ones(1,3),rand(4998,1)+0.2).^5)*10-12)'; // (fictious data on 5000 days)
      
      // Generation of ETP daily time series (5000 days), latitude = 45deg, start = 1/1/1980
      T   = [0.687;0.498;2.774;6.086;10.565;13.702;16.159;15.585;12.619;8.486;3.300;0.778];
      ETP = c_ETP(5000,24,%pi/4,T,198001010000,1);
      
      // GR4J parameters
      X   =[665;1.18;90;3.8;0];
      
      // GR4J simulation
      [Qsim,Ssim,Rsim]=c_GR4J(24,X,P,ETP,[0.6;0.7]);
      
      // Plots
      x=(1:size(P,1))';
      subplot(3,1,1),plot(x,Qsim); // Calculated discharge
      subplot(3,1,2),plot(x,Ssim); // Production store filling level
      subplot(3,1,3),plot(x,Rsim); // Routing store filling level    

This produces the following figure.

HYDROGR

25th of September 2010: Financial Module

The module is dedicated to finance. There are three main areas that are covered: (i) risk measure and management, (ii) asset allocation, and (iii) pricing.

The author of this module is Francesco Menoncin from Brescia University - Economics Department.

For what concerns the risk measure, some functions are dedicated to the computation of Value at Risk (VaR) and Expected Shortfall (ES). Backtest is also implemented in order to check the goodness of such risk measures. Both VaR and ES are also computed in an Extreme Value Theory framework (EVT). Furthermore, it is possible to estimate the parameters of the EVT density function (through maximum likelihood). The Mean Excess Function for graphical study of an EVT distribution is also implemented. The interest rate risk is faced by functions aimed at computing duration, convexity, and yield to maturity. Furthermore, Merton, Vasicek and Cox, Ingersoll and Ross interest rate models are implemented together with the estimation of their parameters. Parametric interpolation of the interest rate curve is possible through both Svennson’s and Nelson-Siegel’s models. Finally, some technical analysis indicators are implemented: Bollinger bands, moving averages, Hurst index.

The asset allocation problem is faced by two functions which compute: (i) the optimal portfolio minimizing the variance of its return and (ii) the optimal portfolio minimizing the expected shortfall of its return. In both cases, the portfolios with and without a riskless asset and with and without short selling are computed.

Pricing problem is approached through functions aimed at: (i) computing the spread on Interest Rate Swaps, (ii) computing the value of options in the Black and Scholes framework (with Greeks and implied volatility), (iii) simulating stochastic processes (through Euler discretization).

Features

The following is an example of the bsgreeks function.

We compute the Greeks on both a call and put option with: underlying price 25 euros, strike price 25 euros, 0.001 (annual) riskless interest rate, 3 month time to maturity (i.e. T=3/12), and 0.2 (annual) volatility.

-->[D,G,Th,R,V]=bsgreeks(25,25,0.01,3/12,0.2)
V =
 4.9727729
R =
 3.0550246
- 3.1793699
Th =
 2.1113101
1.8619344
G =
 4.9727729
D =
 0.5298926
- 0.4701074

public: New Scientific Features in 2010 (last edited 2012-07-16 13:08:34 by michael.baudin@contrib.scilab.org)