function x=dspace(p,sampsiz) ;
% this function simulates a discrete state space process
                                % p is the transition matrix
                                % sampsiz is the length of simulation

[n,n0]=size(p) ;                % determines the number of states = n

% compute the unconditional probability vector which must satisfy
% the condition that pu=p'pu where p is the transition matrix and pu
% is the unconditional probability vector.  This is achieved by
% starting at a guess pu(0) and creating pu(1)=p'pu(0).  The process is
% continued until pu(k) is close to pu(k-1)

pu=ones(n,1)./n ;               % initialize unconditional probabilities
                                % to = 1/n for all j

md=1 ;                          % initialize a convergence criterion
while md > 1e-12 ;              % loop while criterion is large
 pup=p'*pu ;
 md=max(abs(pup-pu)) ;          % compute maximum change in pu
 pu=pup ;                       % update
end ;

ps=cumsum(p')' ;                % computes a matrix ps, whose ji element
                                % is Prob[ y <= y(i) | x=x(j) ] by
                                % summing Prob[ y=y(k) | x=x(j) ]  for
                                % y(k) <= y(i)

pus=cumsum(pu) ;                % computes a vector pus, whose j element
                                % is Prob[ y <= y(j) ] by summing
                                % Prob[ y=y(k) ]  for y(k) <= y(j)

rand('uniform') ;               % initialize random number generator so
                                % that rand returns a number from the
                                % uniform distribution on [0,1]

ib=zeros(2*sampsiz,1) ;         % initializes a vector whose elements
                                % will represent the state index j at
                                % time t - i.e. if we are in the 3rd
                                % state at time t, then ib(t)=3

% the random index generator works as follows.  On the interval [0,1]
% mark the probabilities pus(j) on the axis as in the diagram below for
% the case of n=5.  Draw a uniform random number using the rand function
% in MATLAB.  Suppose it lands between pus(2) and pus(1) as indicated
% below.  Then it makes sense to allocate any draw like that as a case
% in which state 2 is realized since the probability of state 2 is
% pu(2)=pus(2)-pus(1)
%                 rand
%                   |
%                   |
%                   V
% |----------|-------------|---------------|-------------|----------|
% 0        pus(1)        pus(2)          pus(3)        pus(4)   pus(5)=1

% in the next line the find(pus > rand) statement creates a vector whose
% elements are the row numbers j from the vector pus for which pus(j) >
% rand.  So for the example above this would be a column vector with
% elements 2,3,4 and 5.  It then finds the minimum index from among
% these, which in our example is 2, and sets that as the state at time 1

ib(1,1)=min(find(pus > rand)) ; % get initial random index.

% now we loop over the remaining time periods using the conditional
% probabilities instead of the unconditional probabilities.  The
% relevant conditional probability vector at time i, is the row from
% p which corresponds to the state at time i-1

for i=2:(2*sampsiz) ;
 ib(i,1)=min(find(ps(ib(i-1,1),:)' > rand)) ;
end ;

x=ib(sampsiz+1:2*sampsiz,1) ;   % extract the last half of the sample
