SMCSIM Simulates a Switching Markov chain with state space.
Calculates the stationary distribution for a Markov chain. | |
Simulates a Markov chain. | |
Display message and abort function. | |
Kronecker tensor product. | |
Display warning message; disable or enable warning messages. |
Quick test of the routines in module 'markov' |
001 function [x,z] = simsmc(P,Qc,T,init); 002 003 004 005 %SMCSIM Simulates a Switching Markov chain with state space. 006 007 % 008 009 % CALL: [x,z] = smcsim(P,Q,T); 010 011 % [x,z] = smcsim(P,Q,T,init); 012 013 % 014 015 % x = Simulated switching Markov chain 016 017 % z = Simulated Regime process 018 019 % 020 021 % P = Transition matrix for regime process [rxr] 022 023 % Q = Cell array of transition matrices {r,1} 024 025 % Q{i} = Transition matrix for Markov chain i [nxn] 026 027 % T = Length of simulation. 028 029 % init.x0 = Initial state of process x. If not given, it will start from 030 031 % the stationary distribution of minima given z(1). 032 033 % init.z0 = Initial state of regime process. If not given, it will start 034 035 % from the stationary distribution of the Markov chain. 036 037 % 038 039 % Simulates a Switching Markov chain with state space {1,2,...,n}. 040 041 % The regime process has state space {1,2,...,r}. 042 043 % 044 045 % Example: Simulation of a switching Markov chain with two regime states. 046 047 % P=[0.98 0.02; 0.05 0.95]; n=16; 048 049 % Q{1} = rand(n,n)*diag(exp(5*((n:-1:1)-1)/n)); 050 051 % Q{2} = rand(n,n)*diag(exp(5*((1:n)-1)/n)); % They will be normalized to row sum 1. 052 053 % [x,z] = smcsim(P,Q,400); hmmplot(x,z) 054 055 % init.z0 = 2; init.x0 = []; 056 057 % [x,z] = smcsim(P,Q,400,init); hmmplot(x,z,[],[1 2],'','',1) 058 059 % init.z0 = []; init.x0 = 6; 060 061 % [x,z] = smcsim(P,Q,400,init); hmmplot(x,z,[],[1 2],'','',1) 062 063 % Example: Simulation of a Markov chain 064 065 % P=[0.9 0.1; 0.05 0.95]; 066 067 % x = smcsim(1,P,1000); 068 069 % plot(x) 070 071 072 073 % Tested on Matlab 5.3 074 075 % 076 077 % History: 078 079 % Revised by PJ 19-May-2000 080 081 % updated for WAFO 082 083 % Corrected method for simulating starting conditions. 084 085 % Created by PJ (Pär Johannesson) 1997 086 087 % Copyright (c) 1997 by Pär Johannesson 088 089 % Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997 090 091 092 093 % Check input arguments 094 095 096 097 ni = nargin; 098 099 no = nargout; 100 101 error(nargchk(3,4,ni)); 102 103 104 105 if ni < 4, init = []; end 106 107 108 109 if isempty(init) 110 111 init.x0 = []; 112 113 init.z0 = []; 114 115 end 116 117 118 119 % Set constants 120 121 Zstr = '123456789'; 122 123 124 125 126 127 r = length(P); % Number of regime states 128 129 n = length(Qc{1}); % Number of levels 130 131 132 133 % Check that the rowsums of P are equal to 1 134 135 136 137 sumP = sum(P'); 138 139 if sum(sumP == 1) ~= length(P) 140 141 warning(['Rowsums of P not equal to 1. Renormalizing.']); 142 143 for i = 1:length(P) 144 145 P(i,:) = P(i,:)/sumP(i); 146 147 end 148 149 end 150 151 152 153 % Check that the rowsums of Qc{1},...,Qc{r} are equal to 1 154 155 156 157 for i = 1:r 158 159 sumQi = sum(Qc{i}'); 160 161 if sum(sumQi == 1) ~= length(Qc{i}) 162 163 warning(['Rowsums of Q{' Zstr(i) '} not equal to 1. Renormalizing.']); 164 165 for j = 1:length(Qc{i}) 166 167 Qc{i}(j,:) = Qc{i}(j,:)/sumQi(j); 168 169 end 170 171 end 172 173 end 174 175 176 177 178 179 % Make the transition matrix Q for the joint MC (X_k,Z_k) 180 181 182 183 Q = zeros(n*r,n*r); 184 185 I = 0:r:(n-1)*r; 186 187 for z=1:r 188 189 QQ = kron(Qc{z},P); 190 191 Q(I+z,:) = QQ(I+z,:); 192 193 end 194 195 196 197 198 199 % Stationary distribution = ro of Q 200 201 202 203 ro = mc2stat(Q); 204 205 206 207 % Generate random numbers 208 209 210 211 e=rand(T,1); 212 213 214 215 % Start values 216 217 218 219 e0 = e(1); 220 221 if isempty(init.z0) & isempty(init.x0) 222 223 x0z0 = min(find( e0<=cumsum(ro) )); 224 225 x0 = floor((x0z0+1)/r); 226 227 z0 = mod(x0z0-1,r)+1; 228 229 elseif isempty(init.x0) 230 231 z0 = init.z0; 232 233 rox0 = ro(z0:r:end); % Pick stat. distr. for regime z0 234 235 rox0 = rox0/sum(rox0); 236 237 x0 = min(find( e0<=cumsum(rox0) )); 238 239 elseif isempty(init.z0) 240 241 x0 = init.x0; 242 243 z0 = []; % Start from stat. distr of P 244 245 else % Both z0 znd x0 are given 246 247 x0 = init.x0; 248 249 z0 = init.z0; 250 251 end 252 253 254 255 % Simulate Regime process 256 257 258 259 z = mcsim(P,T,z0); 260 261 262 263 % Simulate switching Markov chain 264 265 266 267 x=zeros(T,1); 268 269 x(1) = x0; % First value 270 271 272 273 for k=2:T 274 275 Pi = Qc{z(k)}; 276 277 % eval(['Pi = P' Zstr(z(k)) ';']); 278 279 cumsumPi = cumsum(Pi')'; 280 281 x(k) = min(find( e(k)<=cumsumPi(x(k-1),:) )); 282 283 end 284 285
Comments or corrections to the WAFO group