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";
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]]);
> generator(P);
Warning, No completely valid generator found
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);
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
An example from the Israel, Rosenthal and Wei paper where the matrix calculated using the principal branch of 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]]);
> 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 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]]);
> 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
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]]);
> 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]
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]]);
> generator(P);
generator: Searching for generator with Digits=14
generator/generator: Determinant exceeds product of diagonal elements
See also:
exponential , linalg , LinearAlgebra , Matrix