Download:

Function: generator - find generator of a continuous-time Markov chain with a given time-1 transition matrix

Calling sequence:

generator(P,digits);

Parameters:

P - the transition matrix: a square matrix or Matrix with nonnegative entries and row-sums 1

digits - (optional) positive integer

Description:

Examples:

The function generator must be defined by reading in the file gen6.txt .

> read "gen6.txt";

generator := proc (Po) local P, N, odigs; option `C...

`generator/generator` := proc (Po, Di, odi) local P...

`generator/nextk` := proc (k::evaln, Kmin, Kmax) lo...

`generator/nextki` := proc (k) option `Copyright (c...

The example of Jarrow et al.

> P:= Matrix([[0.8910, 0.0963, 0.0078, 0.0019, 0.0030, 0.0000, 0.0000, 0.0000],
[0.0086, 0.9010, 0.0747, 0.0099, 0.0029, 0.0029, 0.0000, 0.0000],
[0.0009, 0.0291, 0.8894, 0.0649, 0.0101, 0.0045, 0.0000, 0.0009],
[0.0006, 0.0043, 0.0656, 0.8427, 0.0644, 0.0160, 0.0018, 0.0045],
[0.0004, 0.0022, 0.0079, 0.0719, 0.7764, 0.1043, 0.0127, 0.0241],
[0.0000, 0.0019, 0.0031, 0.0066, 0.0517, 0.8246, 0.0435, 0.0685],
[0.0000, 0.0000, 0.0116, 0.0116, 0.0203, 0.0754, 0.6493, 0.2319],
[0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 1.0000]]);

P := _rtable[136858896]

> generator(P);

Warning, No completely valid generator found

_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]
_rtable[138781964]

Note the small negative off-diagonal entries, which are the cause of the warning message. The matrix is a bit easier to see rounded to 4 digits. We also increase infolevel[generator] to see some information about the process.

> infolevel[generator]:= 3;
generator(P,4);

infolevel[generator] := 3

generator:   Searching for generator with Digits=14

generator/generator:   Eigenvalues    [.918611552792567120, .981328562499740320, .884529614120531438, .855134608958769138, .798512366648125526, .704578752690526500, .632061928804652062, 1.]

generator/generator:   k=   []

generator/generator:   Approximate generator with qmin =   -.41983177641404e-3

Warning, No completely valid generator found

_rtable[138790444]

An example from the Israel, Rosenthal and Wei paper where the matrix calculated using the principal branch of ln has some negative off-diagonal elements. In this case the matrix does not have distinct eigenvalues, and generator does not find a valid generator (even though one exists).

> P:= Matrix([
[.284779445, .284035268, .283826586, .147358701],
[.284191780, .284779445, .284035268, .146993507],
[.283487477, .284191780, .284779445, .147541298],
[.284543931, .283487477, .284191780, .147776812]]);

P := _rtable[138795660]

> generator(P);

generator:   Searching for generator with Digits=14

generator/generator:   Eigenvalues    [.99999999999999956, .746358024772408870e-6, .105720032098761312e-2+.181023128288310398e-9*I, .105720032098761312e-2-.181023128288310398e-9*I]

Warning, Matrix does not have distinct eigenvalues

generator/generator:   Going to   25   digits

generator/generator:   Eigenvalues    [.99999999999999956, .746358024772408870e-6, .105720032098761312e-2+.181023128288310398e-9*I, .105720032098761312e-2-.181023128288310398e-9*I]

Warning, Matrix does not have distinct eigenvalues

generator/generator:   k=   [0]

generator/generator:   k=   [-1]

generator/generator:   k=   [1]

generator/generator:   k=   [-2]

generator/generator:   k=   [2]

generator/generator:   k=   [-3]

generator/generator:   k=   [3]

generator/generator:   k=   [-4]

generator/generator:   k=   [4]

Warning, No completely valid generator found

Another example from Israel, Rosenthal and Wu, having more than one valid generator. Since this example also does not have distinct eigenvalues, generator only finds one valid generator (which is Q[2] of the paper). We use the environment variable EnvAllGenerators to search for more than one generator.

> EnvAllGenerators:= true;
b:= exp(-4*Pi);
P:= Matrix([[(2+3*b)/5, (2-2*b)/5, (1-b)/5],
[(2-2*b)/5, (2+3*b)/5, (1-b)/5],
[(2-2*b)/5, (2-2*b)/5, (1+4*b)/5]]);

EnvAllGenerators := true

b := exp(-4*Pi)

P := _rtable[138797956]

> generator(P);

generator:   Searching for generator with Digits=14

generator/generator:   Eigenvalues    [.348734234995040638e-5, 1., .348734234997816196e-5]

Warning, Matrix does not have distinct eigenvalues

generator/generator:   Going to   24   digits

generator/generator:   Eigenvalues    [.348734234995040638e-5, 1., .348734234997816196e-5]

Warning, Matrix does not have distinct eigenvalues

generator/generator:   k=   []

generator/generator:   Possible generator with qmin =   2.51327412322885468

_rtable[138800332]

Here is an example where the eigenvalues are distinct, and generator finds four valid generators.

> P:= Matrix([
[.3333333237, .3333333354, .3333333409],
[.3333333409, .3333333237, .3333333354],
[.3333333354, .3333333409, .3333333237]]);

P := _rtable[138804748]

> generator(P);

generator:   Searching for generator with Digits=14

generator/generator:   Eigenvalues    [.99999999999999988, -.144499999589931634e-7+.476313971357199988e-8*I, -.144499999589931634e-7-.476313971357199988e-8*I]

generator/generator:   Going to   26   digits

generator/generator:   Eigenvalues    [.99999999999999988, -.144499999589931634e-7+.476313971357199988e-8*I, -.144499999589931634e-7-.476313971357199988e-8*I]

generator/generator:   Going to   28   digits

generator/generator:   Eigenvalues    [.99999999999999988, -.144499999589931634e-7+.476313971357199988e-8*I, -.144499999589931634e-7-.476313971357199988e-8*I]

generator/generator:   k=   [0]

generator/generator:   Possible generator with qmin =   4.37036894398100272

generator/generator:   k=   [-1]

generator/generator:   Possible generator with qmin =   4.00269728699146121

generator/generator:   k=   [1]

generator/generator:   Possible generator with qmin =   .742770163046048371

generator/generator:   k=   [-2]

generator/generator:   Possible generator with qmin =   .375098525705631846

generator/generator:   k=   [2]

generator/generator:   k=   [-3]

generator/generator:   k=   [3]

generator/generator:   k=   [-4]

generator/generator:   k=   [4]

generator/generator:   k=   [-5]

generator/generator:   k=   [5]

generator/generator:   k=   [-6]

_rtable[138805412], _rtable[136774684], _rtable[138...
_rtable[138805412], _rtable[136774684], _rtable[138...

An example where the matrix violates one of the necessary conditions for a generator to exist.

> P:= Matrix([[.1,.8,.1],[.1,.1,.8],[.8,.1,.1]]);

P := _rtable[138817156]

> generator(P);

generator:   Searching for generator with Digits=14

generator/generator:   Determinant exceeds product of diagonal elements

See also:

exponential , linalg , LinearAlgebra , Matrix