function [y,BETA,SIGMA] = VARsim()

load BETA.mat;
load SIGMA.mat;

if nargin==0
    T = 190;
    M = 3;
    p = 1;
    BETA = BETA(:,1:M+M*M);
    K = size(BETA,2);
    
    y = zeros(T+p+50,M);
    y(1:p,:) = rand(p,M);
    
    for t = p+1:T+50        
        if t<=50
            BETA(1,:) = BETA(1,:) + randn(1,K)*(0.02*eye(K));
            SIGMA(1,:,:) = squeeze(SIGMA(1,:,:)) + randn(M,M)*(0.01*eye(M));
            bb1 = BETA(1,1:M);
            bb2 = BETA(1,M+1:end)';
            splace = 0; bb3 = 0;
            for ii = 1:p                                          
                for iii = 1:M           
                    bb3(iii,(ii-1)*M+1:ii*M) = bb2(splace+1:splace+M,1)';
                    splace = splace + M;
                end
            end            
            y(t,:) = y(t-1,:)*bb3' + randn(1,M)*squeeze(SIGMA(1,:,:))';
        elseif t>50
            BETA(t-50,:) = BETA(t-50,:) + randn(1,K)*(0.02*eye(K));
            SIGMA(t-50,:,:) = squeeze(SIGMA(t-50,:,:)) + randn(M,M)*(0.01*eye(M));
            bb1 = BETA(t-50,1:M);
            bb2 = BETA(t-50,M+1:end)';
            splace = 0; bb3 = 0;
            for ii = 1:p                                          
                for iii = 1:M           
                    bb3(iii,(ii-1)*M+1:ii*M) = bb2(splace+1:splace+M,1)';
                    splace = splace + M;
                end
            end
            y(t,:) = y(t-1,:)*bb3' + randn(1,M)*squeeze(SIGMA(t-50,:,:));
        end
    end
    BETA = BETA(:,M+1:end);
    y = y(50+1:end,:);
end