clear; close; clc;
% mypath = '/home/user/l571300/private/Ricerca/PhD_batch/large_TVP_MCARLO/';
% mypath = '\\osiride-fs\l571300\private\Ricerca\PhD_batch\large_TVP_MCARLO\';
% addpath(genpath('/Users/fabrizio/Dropbox/PhD/large_TVP/'));
curpath = pwd; % '\\gimecb01\data\ECB business areas\DGI\Databases and Programme files\IPA\Personal folders\Venditti\Research\VARCodes\LargeTVPKMV\';
addpath(genpath(curpath));
mypath = [curpath,'\Out_Mcarlo\'];
mkdir(mypath);

common_setup;
theta_l =  0;
ar(2)   =  0;
prob_break    = 0.025*10; % probability of a break, this implies 1 every year on average, 10 times as high as DGP2.
    for obs = 200  % number of time periods in the exercise
        for m = [7 15]
            for hlow   = theta_l % lower bound for the eigenvalues of the VAR leading matrix. They will fluctuate between this and 0.99
                % some preliminary computations
                mcarlo_prelim;
                for kk=1:Nsimul
                    disp(['m=',num2str(m),' simul = ',num2str(kk)]);                    
                    tic
                    % SIMULATE TVAR COEFFFICIENTS, THIS DIFFERS ACROSS DGPs
                    e=recserar(randn(obs,ncoefs),zeros(1,ncoefs),ar(1)*ones(1,ncoefs)); % shocks for stoch coefficients
                    % plot(e);
                    xbeta=recserar(e,rand(1,ncoefs),ones(1,ncoefs)); % random walk coefficients
                    aa = repmat((maxc(abs(xbeta))),1,obs)';
                    xbeta=(1./aa).*xbeta;                            % bound the RW coefficients to 1
                    
                    % common block of simulations for all DGPs (errors and
                    % volatilities)
                    common_simul;
                    
                    j = 2;
                    beta_shad = NaN*ones(obs,m^2);
                    beta      = beta_shad;                   
                    while j<=obs+1;
                        [betaa,P,L]=matrixaa_rwalkeig(xbeta(j-1,:),theta_l);
                        %% allow for a change in coefficients with probability prob_break
                        beta_shad(j-1,:)=reshape(betaa,1,m^2);                    
                        if j==2||rand<(prob_break)
                            beta(j-1,:) = beta_shad(j-1,:);
                        else
                            beta(j-1,:) = beta(j-2,:);
                        end
                        beta_mat = reshape(beta(j-1,:),m,m);
%                        firstlag(j-1,:)=diag(betaa);
                        y1(j,:)=const*cbeta(j-1,:)+(beta_mat*y1(j-1,:)')'+u(j,:);
                        j=j+1;
                    end;
                    ydata = normaliz(y1(2:end,:)); %% NORMALIZE AND THROW AWAY INITIAL CONDITION
                    kernshape     = {'gauss'};
                    
                    % here estimate models
                    common_estimate;
                    toc;
                    
                end
                if ar(2)==0
                    filename = [mypath,'DGP2bis_m_',num2str(m),'_T_',num2str(obs),'_thetalow_',num2str(theta_l),'.mat'];
                else
                    filename = [mypath,'DGP2bisAR_m_',num2str(m),'_T_',num2str(obs),'_thetalow_',num2str(theta_l),'.mat'];
                end
                save(filename,'kern_dm*','KK_dm*');
                disp(filename)
                clear kern_dm* KK_dm*
            end     % end obs  = obs_grid
        end         % end hlow = theta_l
    end             % end m    = m_grid