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) See also mc2stat, smcsim
Calculates the stationary distribution for a Markov chain. |
Simulates a switching ARMA-process. | |
Simulates a Switching Markov chain with state space. | |
Simulates a switching Markov chain of turning points, | |
Quick test of the routines in module 'markov' |
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 % 021 % See also mc2stat, smcsim 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
Comments or corrections to the WAFO group