Home > wafo > wsim > mcsim.m

# mcsim

## PURPOSE

Simulates a Markov chain.

x=mcsim(P,T,x0)

## DESCRIPTION

``` MCSIM   Simulates a Markov chain.

CALL: x = mcsim(P,T)
x = mcsim(P,T,x0)

x  = Simulated Markov chain.

P  = Transition matrix.    [rxr]
T  = Length of simulation.
x0 = Initial state. (Default:  x0 will be from the stationary
distribution of the Markov chain.)

Simulates a Markov chain with state space {1,2,...,r}

Example:
P=[0.8 0.2; 0.1 0.9];
x=mcsim(P,100); plot(x)

## CROSS-REFERENCE INFORMATION

This function calls:
 mc2stat Calculates the stationary distribution for a Markov chain.
This function is called by:
 sarmasim Simulates a switching ARMA-process. smcsim Simulates a Switching Markov chain with state space. smctpsim Simulates a switching Markov chain of turning points, test_markov Quick test of the routines in module 'markov'

## SOURCE CODE

```001 function x=mcsim(P,T,x0)
002
003 %MCSIM   Simulates a Markov chain.
004 %
005 % CALL: x = mcsim(P,T)
006 %       x = mcsim(P,T,x0)
007 %
008 % x  = Simulated Markov chain.
009 %
010 % P  = Transition matrix.    [rxr]
011 % T  = Length of simulation.
012 % x0 = Initial state. (Default:  x0 will be from the stationary
013 %      distribution of the Markov chain.)
014 %
015 % Simulates a Markov chain with state space {1,2,...,r}
016 %
017 % Example:
018 %   P=[0.8 0.2; 0.1 0.9];
019 %   x=mcsim(P,100); plot(x)
020 %
022
023 % Tested  on Matlab  5.3
024 %
025 % History:
026 % Revised by PJ 26-Jul-2000
027 %   updated help text.
028 % Created by PJ (Pär Johannesson) 1997
029 %   Copyright (c) 1997 by Pär Johannesson
030 %   Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997
031
032 if nargin < 3
033   x0 = [];
034 end
035
036 % Check that the rowsums of P are equal to 1
037
038 sumP = sum(P');
039 if sum(sumP == 1) ~= length(P)
040   disp(['Warning: Rowsums of P not equal to 1. Renormalizing.']);
041   for i = 1:length(P)
042     P(i,:) = P(i,:)/sumP(i);
043   end
044 end
045
046 x=ones(T,1);
047 e=rand(T,1);
048 cumsumP = cumsum(P')';
049
050 % Initial state
051
052 if isempty(x0)    % From stationary distribution
053
054   ro=mc2stat(P);
055
056   x(1) = min(find( e(1)<=cumsum(ro) ));
057 %  x(1) = sum( cumsum(ro) <= e(1) ) + 1;
058
059 else
060
061   x(1) = x0;   % Given
062
063 end
064
065 %---------------------
066
067 for k=2:T
068
069   x(k) = min(find( e(k)<=cumsumP(x(k-1),:) ));
070
071 end
072
073 %  x(k) = sum( cumsumP(x(k-1),:) <= e(k) ) + 1;
074```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003