%______________________ The main program to reproduce Figures 1-6 and Tables 3-5 in Casarin, Grassi, Ravazzolo and Van Dijk (2013)___________________________
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Reproduce Figures 1 and 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
addpath('Figure1');
n = 1500;   %%% Number of monte carlo iterations 
iIter = 1000; %%% iteration 
for fixgenrator=0:1
    run MonteCarloIntegration.m
end
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of Figures 1 and 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Reproduce Table 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Table 2: The four simulation exercises are divided in Complete Model Set and InComplete Model Set, and CPU and GPU computations. Results refer to k=1 Complete Model Set and to k=2 Incomplete Model Set following the panels in table.
%%%%%% COMPLETE MODEL SET
%%%%%%%%%%%%%%%%%%%%% CPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
finaltime = NaN(2,1);
addpath('../DeCo_Package/codeCPU');
for k=1:2
    bol.bias=k;    
    % _______________________Simulate true model and predictive ________
    if bol.bias==1
        load('input/BiasedPred.mat');% / is for Linux and Windws compatibility. If is not working under Windows substitute / with \
        dire='output/Biased/';
        vY=data.Y;
        mX=permute(data.fcst,[1 3 2]);
    else
        load('input/UnBiasedPred.mat');% / is for Linux and Windws compatibility. If is not working under Windows substitute / with \
        dire='output/UnBiased/';
        vY=data.Y;
        mX=permute(data.fcst,[1 3 2]);
    end
    % _____________________ Filtering __________________________________
    iL=1;
    [cT, iM, iKL] = size(mX);
    iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
    vY = vY(end - cT + 1:end, iSelectSeries);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                         Parameter setting                               %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Setting.cN = 1000; %5000                            %%% Number of particles
    Setting.iDimOmega = iKL;                            %%% State space dimension
    Setting.iM = iM;                                    %%% M parameter
    Setting.iL = iL;                                    %%% L parameter
    Setting.cT = cT;                                    %%% Time dimension
    Setting.dLambda = 0.95;                             %%% Lambda matrix
    Setting.iEstimate = 0;                           %%% Estimate or not the Lambda and Sigma
    Setting.dKappa = 1;                                 %%% Resampling parameter
    iLearning=0;
    iTau=9;
    if iLearning == 0                                   %%% Learning or not
        Setting.iTau = -1;
    else
        Setting.iTau = iTau;
    end
    %%-------------------- Setting the Sigma and Lambda matrix --------------%%
    mDataRand = [0.01,0.1];
    mSigmaTemp = [];
    vSigmaTemp = [];
    for i = 1:Setting.iL
        vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
        mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
    end
    Setting.Sigma = mSigmaTemp;
    Setting.vSigma = vSigmaTemp;
    %%% ------------ Create the matrix that will save the results ----------%%%
    mWeights = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
    mYPred = zeros(Setting.cT, 3,  Setting.iL, Setting.iM);
    TrueweightsCumRange = zeros(Setting.cT, Setting.iDimOmega, 4);
    mSigma = zeros(Setting.cT, Setting.iL, Setting.iM);
    mLambda = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
    %%% ------------------ Start run filtering and timing ----------------- %%%
    tic;
    parfor g1 = 1:Setting.iM
        mError = zeros(Setting.cT, Setting.iDimOmega);
        mXInit = zeros(Setting.cT, Setting.iDimOmega);
        disp(char('We are at the iteration: ', num2str(g1)));
        %%% ---------------- Construct the mX and mError ------------------ %%%
        mXInit(:, :) = reshape(mX(:, g1,  :), Setting.cT, Setting.iDimOmega);
        mError(:, :) = (Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, 1)) - mXInit(:, :)).^2;
        [S0] = initPSSimCPU(Setting, vY, mXInit);                             %% COMMENT: the initialization in BCRVD(2013) is done using values larger than the priors (and not equal to priors as in the empirical application and in the DeCo toolbox). To replicate BCRVD(2013) results we use the same values 
        [SSta0] = initStaCPU(Setting);
        [S, SSta] = PFCPU(Setting, S0, SSta0, vY, mXInit, mError);
        mWeightsTemp = zeros(Setting.cT, Setting.iDimOmega);
        for i = 1:Setting.cT
            mWeightsTemp(i, :) = logmul(SSta.mXFilt(i, :));
        end
        mWeights(:, :, g1) = mWeightsTemp;
        mYPred(:, :, :, g1) = SSta.mYPred;
        mSigma(:, :, g1) = SSta.mSigmaFilt;
        mLambda(:, :, :, g1) = SSta.mLambdaFilt;
    end
    Time = toc;
    finaltime(k,1)=Time;
    % _____________________ statistics __________________________________
    for kk=1:Setting.iDimOmega
        for tt=1:Setting.cT
            TrueweightsCumRange(tt, kk, 1) = quantile(mWeights(tt, kk, :),0.025);
            TrueweightsCumRange(tt, kk, 2) = median(mWeights(tt, kk, :),3);
            TrueweightsCumRange(tt, kk, 3) = quantile(mWeights(tt, kk, :),0.975);
            TrueweightsCumRange(tt, kk, 4) = mean(mWeights(tt, kk, :), 3);
        end
    end
    save(strcat(dire,'weightsEx1CPU.mat'),'TrueweightsCumRange');
    save(strcat(dire,'PredEx1CPU.mat'),'mYPred');
end
save('output/TimeExercisesTable2CompleteModelSetCPU.mat', 'finaltime');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% GPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
finaltime = NaN(2,1);
addpath('../DeCo_Package/codeGPU');
for kex=1:2
    bol.bias=kex;    
    % _______________________Simulate true model and predictive ________
    if bol.bias==1
        load('input/BiasedPred.mat');
        dire = 'output/Biased/';
        vY = data.Y;
        mX = permute(data.fcst,[1 3 2]);
    else
        load('input/UnBiasedPred.mat');
        dire = 'output/UnBiased/';
        vY = data.Y;
        mX = permute(data.fcst,[1 3 2]);
    end
    iL=1;
    [cT, iM, iKL] = size(mX);
    iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
    vY = vY(end - cT + 1:end, iSelectSeries);
    %                       Parameter setting                                 %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    iDraws=100;
    Setting.cT = cT;                                    %%% Time dimension
    Setting.cN = 1000; %5000                            %%% Number of particles
    Setting.iDimOmega = iKL;                            %%% State space dimension
    Setting.dLambda = 0.95;                             %%% Lambda parameter
    Setting.iL = iL;                                    %%% L parameter
    Setting.iM = iM;                                    %%% M parameter
    Setting.iDraws = iDraws;                               %%% Draws from the predictive
    Setting.dKappa = 1;                                 %%% Resampling parameter
    Setting.iEstimate = 0;                           %%% Estimate or not the Lambda and Sigma
    Setting.dKappa = 1;                                 %%% Resampling parameter
    iLearning=0;
    iTau=9;
    if iLearning == 0                                   %%% Learning or not
        Setting.iTau = -1;
    else
        Setting.iTau = iTau;
    end
    %%-------------------- Setting the Sigma and Lambda matrix --------------%%
    mDataRand=[0.01,0.1];
    mSigmaTemp = [];
    vSigmaTemp = [];
    for i = 1:Setting.iL
        vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
        mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
    end
    Setting.Sigma = mSigmaTemp;
    Setting.vSigma = vSigmaTemp;
    % _____________________ Filtering __________________________________
    weightsCum = zeros(Setting.cT, Setting.iDimOmega, Setting.iDraws);
    TrueweightsCumRange = zeros(Setting.cT, Setting.iDimOmega, 4);
    iCounter = 0;
    vVector = 1:Setting.iDraws:Setting.iM + 1;
    try           %%%--------- Catch possible error on GPU memory ---------%%%
    tic;
    for i=1:size(vVector, 2) - 1
        disp(char('We are at the iteration: ', num2str(i)));
        mSigma = [];
        mXInit = zeros(Setting.cT * Setting.iDraws, Setting.iDimOmega * Setting.iDraws);
        mY = zeros(Setting.cT * Setting.iDraws, Setting.iDraws);
        iTemp = 0; iSeries = 0; iIter = 1;
        for g = (vVector(i):vVector(i + 1) - 1)
            %%% ---------------- Construct the mX and mError -------------- %%%
            mXInit(iTemp + 1:iTemp + Setting.cT, iSeries + 1:iSeries + Setting.iDimOmega) = reshape(mX(:, g, :), Setting.cT, Setting.iDimOmega);
            mY(iTemp + 1:iTemp + Setting.cT, iIter) = vY(:, 1);
            iTemp = iTemp + Setting.cT;
            iSeries = iSeries + Setting.iDimOmega;
            mSigma = blkdiag(mSigma, Setting.Sigma); %% This is the Sigma for the PF
            iIter = iIter + 1;
        end
        %%%%%%%% Load on the GPU %%%%%%%%
        mXGPU = gpuArray(mXInit);
        mYGPU = gpuArray(mY);
        %%%%%%%%%%% Initializzation of the elements %%%%%%%%%%%%%%%%%%%%%%%%%%%
        [S0_GPU] = initPSGPUSim(Setting, mYGPU, mXGPU);                    %% COMMENT: the initialization in BCRVD(2013) is done using values larger than the priors (and not equal to priors as in the empirical application and in the DeCo toolbox). To replicate BCRVD(2013) results we use the same values 
        %%%%%%%%%%%%%%%% Clean the memory before Particle Filter %%%%%%%%%%%%%%
        clear mYGPU mXGPU;
        %%%%%%%%%%%%%%%%%%%%%%% Initialize the states %%%%%%%%%%%%%%%%%%%%%%%%%
        [SSta0_GPU] = initStaGPU(Setting);
        %%%%%%%%%%%%%%%% Preparing dataset for PF   %%%%%%%%%%%%%%%%%%%%%%%%%%%
        mXTot = zeros(Setting.iDimOmega * Setting.iDraws, Setting.cT, Setting.iL);
        mError = zeros(Setting.cT, Setting.iDimOmega * Setting.iDraws, Setting.iL);
        z = 1;
        for k = 1:Setting.iDraws
            mXTot(z:z + Setting.iDimOmega - 1, :) = reshape(mX(:, k + iCounter, :), Setting.cT, Setting.iDimOmega)';   
            vYTemp = kron(ones(1, Setting.iDimOmega), vY(:, 1));    %%%%Reshape the observations to create the learning errors
            vYTemp = Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, 1));    %%%%Reshape the observations to create the learning errors           
            mError(:, z:z + Setting.iDimOmega - 1) =  (vYTemp - reshape(mX(:, k + iCounter, :), Setting.cT, Setting.iDimOmega)).^2; %%%  Construct the errors
            z = z + Setting.iDimOmega;
        end
        %%%%%%%%%%%%%%%%%%%%% Starting the particle filter   %%%%%%%%%%%%%%%%%%
        SettingGPU = struct('cN', gpuArray(Setting.cN),'iDimOmega', gpuArray(Setting.iDimOmega),...
            'Sigma', gpuArray(Setting.Sigma),'vSigma',gpuArray(Setting.vSigma), ...
            'cT', gpuArray(Setting.cT), 'iDraws', gpuArray(Setting.iDraws), ...
            'iTau', gpuArray(Setting.iTau), 'dLambda', gpuArray(Setting.dLambda), ...
            'iL', gpuArray(Setting.iL), 'iM', gpuArray(Setting.iM), 'iEstimate', gpuArray(Setting.iEstimate), 'dKappa', gpuArray(Setting.dKappa));
        %%%%%%%%%%%%%%%%%%%%% Run the Particle filter on the GPU %%%%%%%%%%%%%%
        [~, SSta] = PFGPU(SettingGPU, S0_GPU, SSta0_GPU, gpuArray(vY), gpuArray(mXTot), gpuArray(mError));
        %%%%%%%%%%%%%%%%%% Go back to the CPU and do managing %%$$$$$$%%%%%%%%%
        mXFilt = gather(SSta.mXFilt);
        mYCum = gather(reshape(SSta.mYPred, Setting.cT * Setting.iDraws,  Setting.cN, []));
        clear SSta S S0_GPU SSta0_GPU;     %% Cancel old variables to freee memory in the GPU
        %%%%%%%%%% Managing the mean and quantile of the predictive %%%%%%%%%%%
        for g = 1:Setting.iDraws
            for kk=1:Setting.iL
                mYPred(:, 1, kk, g + iCounter) =  median(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk)), 2);
                mYPred(:, 2, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.025)';
                mYPred(:, 3, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.975)';
            end
        end
        iCounter = iCounter + Setting.iDraws;
        %%%%%%%%%%%%%%%%%%%% Managing the Weights %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for ga = 1:Setting.cT
            z = 1;
            for j = 1:Setting.iDraws
                weightsCum(ga, :, j + vVector(i) - 1) = logmul(squeeze(mXFilt(z:z + Setting.iDimOmega-1, ga))');
                z = z + Setting.iDimOmega;
            end
        end
    end
    Time = toc;
    finaltime(kex,1)=Time;
    % _____________________ statistics __________________________________
    for kk=1:Setting.iDimOmega
        for tt=1:Setting.cT
            TrueweightsCumRange(tt, kk, 1) = quantile(weightsCum(tt, kk, :), 0.025);
            TrueweightsCumRange(tt, kk, 2) = median(weightsCum(tt, kk, :), 3);
            TrueweightsCumRange(tt, kk, 3) = quantile(weightsCum(tt, kk, :), 0.975);
            TrueweightsCumRange(tt, kk, 4) = mean(weightsCum(tt, kk, :), 3);
        end
    end
    catch
       disp('-----------------------------------------------------------------');
       disp(char('The dimension of the particles', num2str(Setting.cN), 'and the dimension of the blocks', num2str(iDraws)));
       disp(char('are too big for the RAM available on the GPU card.')); 
       disp(char('Please try to decrease the dimension of the particles or the blocks or both.'));
       disp('-----------------------------------------------------------------');
   break;
    end
save(strcat(dire,'weightsEx1GPU.mat'),'TrueweightsCumRange');
save(strcat(dire,'PredEx1GPU.mat'),'mYPred');
end
save('output\TimeExercisesTable2CompleteModelSetGPU.mat', 'finaltime');
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% INCOMPLETE MODEL SET
%%%%%%%%%%%%%%%%%%%%%% CPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
finaltime = NaN(2,1);
addpath('../DeCo_Package/codeCPU');
for k=1:2
    bol.bias=k;    
    % _______________________Simulate true model and predictive ________
    if bol.bias==1
        load('input/BiasedPred.mat');
        dire = 'output/Biased/';
        vY = data.Y;
        mX = permute(data.fcst,[1 3 2]);
    else
        load('input/UnBiasedPred.mat');
        dire = 'output/UnBiased/';
        vY = data.Y;
        mX = permute(data.fcst,[1 3 2]);
    end
    mX = mX(:, :, 2:end);
    % _____________________ Filtering __________________________________
    iL=1;
    [cT, iM, iKL] = size(mX);
    iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
    vY = vY(end - cT + 1:end, iSelectSeries);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                         Parameter setting                               %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Setting.cN = 1000; %5000                            %%% Number of particles
    Setting.iDimOmega = iKL;                            %%% State space dimension
    Setting.iM = iM;                                    %%% M parameter
    Setting.iL = iL;                                    %%% L parameter
    Setting.cT = cT;                                    %%% Time dimension
    Setting.dLambda = 0.95;                             %%% Lambda matrix
    Setting.iEstimate = 0;                                %%% Sigma value
    Setting.dKappa = 1;                                 %%% Resampling parameter
    iLearning=0;
    if iLearning == 0                                   %%% Learning or not
        Setting.iTau = -1;
    else
        Setting.iTau = iTau;
    end
    %%-------------------- Setting the Sigma and Lambda matrix --------------%%
    mDataRand=[0.01,0.1];
    mSigmaTemp = [];
    vSigmaTemp = [];
    for i = 1:Setting.iL
        vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
        mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
    end
    Setting.Sigma = mSigmaTemp;
    Setting.vSigma = vSigmaTemp;
    %%% ------------ Create the matrix that will save the results ----------%%%
    mWeights = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
    mYPred = zeros(Setting.cT, 3,  Setting.iL, Setting.iM);
    TrueweightsCumRange = zeros(Setting.cT, Setting.iDimOmega, 4);
    mSigma = zeros(Setting.cT, Setting.iL, Setting.iM);
    mLambda = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
    %%% ------------------ Start run filtering and timing ----------------- %%%
    tic;
    parfor g1 = 1:Setting.iM
        mError = zeros(Setting.cT, Setting.iDimOmega);
        mXInit = zeros(Setting.cT, Setting.iDimOmega);
        disp(char('We are at the iteration: ', num2str(g1)));
        %%% ---------------- Construct the mX and mError ------------------ %%%
        mXInit(:, :) = reshape(mX(:, g1,  :), Setting.cT, Setting.iDimOmega);
        mError(:, :) = (Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, 1)) - mXInit(:, :)).^2;
        [S0] = initPSSimCPU(Setting, vY, mXInit);                             %% COMMENT: the initialization in BCRVD(2013) is done using values larger than the priors (and not equal to priors as in the empirical application and in the DeCo toolbox). To replicate BCRVD(2013) results we use the same values 
        [SSta0] = initStaCPU(Setting);
        [S, SSta] = PFCPU(Setting, S0, SSta0, vY, mXInit, mError);
        mWeightsTemp = zeros(Setting.cT, Setting.iDimOmega);
        for i = 1:Setting.cT
            mWeightsTemp(i, :) = logmul(SSta.mXFilt(i, :));
        end
        mWeights(:, :, g1) = mWeightsTemp;
        mYPred(:, :, :, g1) = SSta.mYPred;
        mSigma(:, :, g1) = SSta.mSigmaFilt;
        mLambda(:, :, :, g1) = SSta.mLambdaFilt;
    end
    Time = toc;
    finaltime(k,1)=Time;
    % _____________________ statistics __________________________________
    for kk=1:Setting.iDimOmega
        for tt=1:Setting.cT
            TrueweightsCumRange(tt, kk, 1) = quantile(mWeights(tt,kk,:),0.025);
            TrueweightsCumRange(tt, kk, 2) = median(mWeights(tt,kk,:),3);
            TrueweightsCumRange(tt, kk, 3) = quantile(mWeights(tt,kk,:),0.975);
            TrueweightsCumRange(tt, kk, 4) = mean(mWeights(tt,kk,:),3);
        end
    end
    save(strcat(dire, 'weightsEx2CPU.mat'), 'TrueweightsCumRange');
    save(strcat(dire, 'PredEx2CPU.mat'), 'mYPred');
end    
save('output\TimeExercisesTable2IncompleteModelSetCPU.mat', 'finaltime');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% GPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
addpath('../DeCo_Package/codeGPU');
finaltime=NaN(2,1);
for kex=1:2
    bol.bias=kex;    
    % _______________________Simulate true model and predictive ________
    if bol.bias==1
        load('input/BiasedPred.mat');
        dire='output/Biased/';
        vY=data.Y;
        mX=permute(data.fcst,[1 3 2]);
    else
        load('input/UnBiasedPred.mat');
        dire='output/UnBiased/';
        vY=data.Y;
        mX=permute(data.fcst,[1 3 2]);
    end
    mX=mX(:,:,2:end);
    % _____________________ Parameter setting __________________________
    iL=1;
    [cT, iM, iKL] = size(mX);
    iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
    vY = vY(end - cT + 1:end, iSelectSeries);
    %                       Parameter setting                                 %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    iDraws=100;
    Setting.cT = cT;                                    %%% Time dimension
    Setting.cN = 1000; %5000                            %%% Number of particles
    Setting.iDimOmega = iKL;                            %%% State space dimension
    Setting.dLambda = 0.95;                             %%% Lambda parameter
    Setting.iL = iL;                                    %%% L parameter
    Setting.iM = iM;                                    %%% M parameter
    Setting.iDraws = iDraws;                               %%% Draws from the predictive
    Setting.dKappa = 1;                                 %%% Resampling parameter
    Setting.iEstimate = 0;                           %%% Estimate or not the Lambda and Sigma
    Setting.dKappa = 1;                                 %%% Resampling parameter
    iLearning=0;
    iTau=9;
    if iLearning == 0                                   %%% Learning or not
        Setting.iTau = -1;
    else
        Setting.iTau = iTau;
    end
    %%-------------------- Setting the Sigma and Lambda matrix --------------%%
    mDataRand=[0.01,0.1];
    mSigmaTemp = [];
    vSigmaTemp = [];
    for i = 1:Setting.iL
        vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
        mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
    end
    Setting.Sigma = mSigmaTemp;
    Setting.vSigma = vSigmaTemp;
    % _____________________ Filtering __________________________________
    weightsCum = zeros(Setting.cT, Setting.iDimOmega, Setting.iDraws);
    TrueweightsCumRange = zeros(Setting.cT, Setting.iDimOmega, 4);
    iCounter = 0;
    vVector = 1:Setting.iDraws:Setting.iM + 1;
    try
    tic;
    for i=1:size(vVector, 2) - 1
        disp(char('We are at the iteration: ', num2str(i)));
        mSigma = [];
        mXInit = zeros(Setting.cT * Setting.iDraws, Setting.iDimOmega * Setting.iDraws);
        mY = zeros(Setting.cT * Setting.iDraws, Setting.iDraws);
        iTemp = 0; iSeries = 0; iIter = 1;
        for g = (vVector(i):vVector(i + 1) - 1)
            %%% ---------------- Construct the mX and mError -------------- %%%
            mXInit(iTemp + 1:iTemp + Setting.cT, iSeries + 1:iSeries + Setting.iDimOmega) = reshape(mX(:, g, :), Setting.cT, Setting.iDimOmega);
            mY(iTemp + 1:iTemp + Setting.cT, iIter) = vY(:, 1);
            iTemp = iTemp + Setting.cT;
            iSeries = iSeries + Setting.iDimOmega;
            mSigma = blkdiag(mSigma, Setting.Sigma); %% This is the Sigma for the PF
            iIter = iIter + 1;
        end
        %%%%%%%% Load on the GPU %%%%%%%%
        mXGPU = gpuArray(mXInit);
        mYGPU = gpuArray(mY);
        %%%%%%%%%%% Initializzation of the elements %%%%%%%%%%%%%%%%%%%%%%%%%%%
        [S0_GPU] = initPSGPUSim(Setting, mYGPU, mXGPU);                    %% COMMENT: the initialization in BCRVD(2013) is done using values larger than the priors (and not equal to priors as in the empirical application and in the DeCo toolbox). To replicate BCRVD(2013) results we use the same values
        %%%%%%%%%%%%%%%% Clean the memory before Particle Filter %%%%%%%%%%%%%%
        clear mYGPU mXGPU;
        %%%%%%%%%%%%%%%%%%%%%%% Initialize the states %%%%%%%%%%%%%%%%%%%%%%%%%
        [SSta0_GPU] = initStaGPU(Setting);
        %%%%%%%%%%%%%%%% Preparing dataset for PF   %%%%%%%%%%%%%%%%%%%%%%%%%%%
        mXTot = zeros(Setting.iDimOmega * Setting.iDraws, Setting.cT, Setting.iL);
        mError = zeros(Setting.cT, Setting.iDimOmega * Setting.iDraws, Setting.iL);
        z = 1;
        for k = 1:Setting.iDraws
            mXTot(z:z + Setting.iDimOmega - 1, :) = reshape(mX(:, k + iCounter, :), Setting.cT, Setting.iDimOmega)';  
            vYTemp = kron(ones(1, Setting.iDimOmega), vY(:, 1));    %%%%Reshape the observations to create the learning errors
            vYTemp = Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, 1));    %%%%Reshape the observations to create the learning errors 
            mError(:, z:z + Setting.iDimOmega - 1) =  (vYTemp - reshape(mX(:, k + iCounter, :), Setting.cT, Setting.iDimOmega)).^2; %%%  Construct the errors
            z = z + Setting.iDimOmega;
        end
        %%%%%%%%%%%%%%%%%%%%% Starting the particle filter   %%%%%%%%%%%%%%%%%%
        SettingGPU = struct('cN', gpuArray(Setting.cN),'iDimOmega', gpuArray(Setting.iDimOmega),...
            'Sigma', gpuArray(Setting.Sigma),'vSigma',gpuArray(Setting.vSigma), ...
            'cT', gpuArray(Setting.cT), 'iDraws', gpuArray(Setting.iDraws), ...
            'iTau', gpuArray(Setting.iTau), 'dLambda', gpuArray(Setting.dLambda), ...
            'iL', gpuArray(Setting.iL), 'iM', gpuArray(Setting.iM), 'iEstimate', gpuArray(Setting.iEstimate), 'dKappa', gpuArray(Setting.dKappa));
        %%%%%%%%%%%%%%%%%%%%% Run the Particle filter on the GPU %%%%%%%%%%%%%%
        [~, SSta] = PFGPU(SettingGPU, S0_GPU, SSta0_GPU, gpuArray(vY), gpuArray(mXTot), gpuArray(mError));
        %%%%%%%%%%%%%%%%%% Go back to the CPU and do managing %%$$$$$$%%%%%%%%%
        mXFilt = gather(SSta.mXFilt);
        mYCum = gather(reshape(SSta.mYPred, Setting.cT * Setting.iDraws,  Setting.cN, []));
        clear SSta S S0_GPU SSta0_GPU;     %% Cancel old variables to freee memory in the GPU
        %%%%%%%%%% Managing the mean and quantile of the predictive %%%%%%%%%%%
        for g = 1:Setting.iDraws
            for kk=1:Setting.iL
                mYPred(:, 1, kk, g + iCounter) =  median(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk)), 2);
                mYPred(:, 2, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.025)';
                mYPred(:, 3, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.975)';
            end
        end
        iCounter = iCounter + Setting.iDraws;
        %%%%%%%%%%%%%%%%%%%% Managing the Weights %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for ga = 1:Setting.cT
            z = 1;
            for j = 1:Setting.iDraws
                weightsCum(ga, :, j + vVector(i) - 1) = logmul(squeeze(mXFilt(z:z + Setting.iDimOmega - 1, ga))');
                z = z + Setting.iDimOmega;
            end
        end
    end
    Time=toc;
    finaltime(kex,1)=Time;
    % _____________________ statistics __________________________________
    for kk=1:Setting.iDimOmega
        for tt=1:Setting.cT
            TrueweightsCumRange(tt, kk, 1) = quantile(weightsCum(tt,kk,:),0.025);
            TrueweightsCumRange(tt, kk, 2) = median(weightsCum(tt,kk,:),3);
            TrueweightsCumRange(tt, kk, 3) = quantile(weightsCum(tt,kk,:),0.975);
            TrueweightsCumRange(tt, kk, 4) = mean(weightsCum(tt,kk,:),3);
        end
    end
    catch
       disp('-----------------------------------------------------------------');
       disp(char('The dimension of the particles', num2str(Setting.cN), 'and the dimension of the blocks', num2str(iDraws)));
       disp(char('are too big for the RAM available on the GPU card.')); 
       disp(char('Please try to decrease the dimension of the particles or the blocks or both.'));
       disp('-----------------------------------------------------------------');
     break;
    end
    save(strcat(dire,'weightsEx2GPU.mat'),'TrueweightsCumRange');
    save(strcat(dire,'PredEx2GPU.mat'),'mYPred');
end
save('output/TimeExercisesTable2IncompleteModelSetGPU.mat', 'finaltime');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of Table 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Reproduce Figure 3-4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
for k=1:2
    bol.bias=k;
    if bol.bias==1
        dire='output/Biased/';
    else
        dire='output/UnBiased/';
    end
    load(strcat(dire,'weightsEx1GPU.mat'));
    GPUweights1=TrueweightsCumRange;
    load(strcat(dire,'weightsEx1CPU.mat'));
    CPUweights1=TrueweightsCumRange;
    
    %% Weights
    figure
    subplot(3, 1, 1)
    plot(squeeze(CPUweights1(:, 1, 2)), '-r');hold on;
    plot(squeeze(GPUweights1(:, 1, 2)), '-b');hold on;
    plot(squeeze(CPUweights1(:, 1, [1,3])), '--r');hold on;
    plot(squeeze(GPUweights1(:, 1, [1,3])), '--b');hold off;
    legend('GPU', 'CPU');
    subplot(3, 1, 2)
    plot(squeeze(CPUweights1(:, 2, 2)), '-r');hold on;
    plot(squeeze(GPUweights1(:, 2, 2)), '-b');hold on;
    plot(squeeze(CPUweights1(:, 2, [1,3])), '--r');hold on;
    plot(squeeze(GPUweights1(:, 2, [1,3])), '--b');hold off;
    legend('GPU', 'CPU');
    subplot(3, 1, 3)
    plot(squeeze(CPUweights1(:, 3, 2)), '-r');hold on;
    plot(squeeze(GPUweights1(:, 3, 2)), '-b');hold on;
    plot(squeeze(CPUweights1(:, 3, [1,3])), '--r');hold on;
    plot(squeeze(GPUweights1(:, 3, [1,3])), '--b');hold off;
    legend('GPU', 'CPU');
    if bol.bias==1
        figname='output/Figure3Biased.pdf';
    else
        figname='output/Figure3UnBiased.pdf';
    end
    saveas(gcf, figname);
    
    %% Predictive density
    load(strcat(dire,'PredEx1GPU.mat'));
    PredGPU1 = mYPred;
    load(strcat(dire,'PredEx1CPU.mat'));
    PredCPU1 = mYPred;
    figure
    mHelpCPU = squeeze(mean(PredCPU1(1:end,:, :), 3));
    mHelpGPU = squeeze(mean(PredGPU1(1:end,:, :), 3));
    plot(mHelpCPU(:, 1), '-r');hold on;
    plot(mHelpGPU(:, 1), '-b');hold on;
    plot(mHelpCPU(:, [2,3]), '--r');hold on;
    plot(mHelpGPU(:, [2,3]), '--b');hold on;
    legend('GPU', 'CPU');
    if bol.bias==1
        figname='output/Figure4Biased.pdf';
    else
        figname='output/Figure4UnBiased.pdf';
    end
    saveas(gcf, figname);
end
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of Figures 3-4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Reproduce Figure 5-6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
for k=1:2
    bol.bias=k;
    if bol.bias==1
        dire='output/Biased/';
    else
        dire='output/UnBiased/';
    end
    load(strcat(dire,'weightsEx2GPU.mat'));
    GPUweights1 = TrueweightsCumRange;
    load(strcat(dire,'weightsEx2CPU.mat'));
    CPUweights1 = TrueweightsCumRange;
    
    %% Weights
    figure
    subplot(2, 1, 1)
    plot(squeeze(CPUweights1(:, 1, 2)), '-r');hold on;
    plot(squeeze(GPUweights1(:, 1, 2)), '-b');hold on;
    plot(squeeze(CPUweights1(:, 1, [1,3])), '--r');hold on;
    plot(squeeze(GPUweights1(:, 1, [1,3])), '--b');hold off;
    legend('GPU', 'CPU');
    subplot(2, 1, 2)
    plot(squeeze(CPUweights1(:, 2, 2)), '-r');hold on;
    plot(squeeze(GPUweights1(:, 2, 2)), '-b');hold on;
    plot(squeeze(CPUweights1(:, 2, [1,3])), '--r');hold on;
    plot(squeeze(GPUweights1(:, 2, [1,3])), '--b');hold off;
    legend('GPU', 'CPU');
    if bol.bias==1
        figname='output/Figure5Biased.pdf';
    else
        figname='output/Figure5UnBiased.pdf';
    end
    saveas(gcf, figname);
    
    %% Predictive density
    load(strcat(dire,'PredEx2GPU.mat'));
    PredGPU1 = mYPred;
    load(strcat(dire,'PredEx2CPU.mat'));
    PredCPU1 = mYPred;
    figure
    mHelpCPU = squeeze(mean(PredCPU1(1:end,:, :), 3));
    mHelpGPU = squeeze(mean(PredGPU1(1:end,:, :), 3));
    plot(mHelpCPU(:, 1), '-r');hold on;
    plot(mHelpGPU(:, 1), '-b');hold on;
    plot(mHelpCPU(:, [2,3]), '--r');hold on;
    plot(mHelpGPU(:, [2,3]), '--b');hold on;
    legend('GPU', 'CPU');
    if bol.bias==1
        figname='output/Figure6Biased.pdf';
    else
        figname='output/Figure6UnBiased.pdf';
    end
    saveas(gcf, figname);
end
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of Figures 5-6 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    Reproduce Table 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% CPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
finaltime = NaN(2,1);
addpath('../DeCo_Package/codeCPU');
for k=1:2
    bol.bias=k;    
    % _______________________Load data and predictive ________
    load('input/UnBiasedPred.mat');
    vY=data.Y;
    if bol.bias==1
        mX=permute(data.fcst,[1 3 2]);
    else
        mX=permute(data.fcst(:,2:end,:),[1 3 2]);
    end
    % _____________________ Filtering __________________________________
    iL=1;
    [cT, iM, iKL] = size(mX);
    iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
    vY = vY(end - cT + 1:end, iSelectSeries);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                         Parameter setting                               %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Setting.cN = 1000; %5000                            %%% Number of particles
    Setting.iDimOmega = iKL;                            %%% State space dimension
    Setting.iM = iM;                                    %%% M parameter
    Setting.iL = iL;                                    %%% L parameter
    Setting.cT = cT;                                     %%% Time dimension
    Setting.dLambda = 0.95;                             %%% Lambda value
    Setting.iEstimate = 0.01;                           %%% Estimate or not the Lambda and Sigma
    Setting.dKappa = 1;                                 %%% Resampling parameter
    iLearning=1;
    iTau=9;
    if iLearning == 0                                   %%% Learning or not
        Setting.iTau = -1;
    else
        Setting.iTau = iTau;
    end
    %%-------------------- Setting the Sigma and Lambda matrix --------------%%
    mDataRand=[0.01,0.1];
    mSigmaTemp = [];
    vSigmaTemp = [];
    for i = 1:Setting.iL
        vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
        mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
    end
    Setting.Sigma = mSigmaTemp;
    Setting.vSigma = vSigmaTemp;
    %%% ------------ Create the matrix that will save the results ----------%%%
    mWeights = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
    mYPred = zeros(Setting.cT, 3,  Setting.iL, Setting.iM);
    mSigma = zeros(Setting.cT, Setting.iL, Setting.iM);
    mLambda = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
    %%% ------------------ Start run filtering and timing ----------------- %%%
    tic;
    parfor g1 = 1:Setting.iM
        mError = zeros(Setting.cT, Setting.iDimOmega);
        mXInit = zeros(Setting.cT, Setting.iDimOmega);
        disp(char('We are at the iteration: ', num2str(g1)));
        %%% ---------------- Construct the mX and mError ------------------ %%%
        mXInit(:, :) = reshape(mX(:, g1,  :), Setting.cT, Setting.iDimOmega);
        mError(:, :) = (Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, 1)) - mXInit(:, :)).^2;
        [S0] = initPSSimCPU(Setting, vY, mXInit);                             %% COMMENT: the initialization in BCRVD(2013) is done using values larger than the priors (and not equal to priors as in the empirical application and in the DeCo toolbox). To replicate BCRVD(2013) results we use the same values
        [SSta0] = initStaCPU(Setting);
        [S, SSta] = PFCPU(Setting, S0, SSta0, vY, mXInit, mError);
        mWeightsTemp = zeros(Setting.cT, Setting.iDimOmega);
        for i = 1:Setting.cT
            mWeightsTemp(i, :) = logmul(SSta.mXFilt(i, :));
        end
        mWeights(:, :, g1) = mWeightsTemp;
        mYPred(:, :, :, g1) = SSta.mYPred;
        mSigma(:, :, g1) = SSta.mSigmaFilt;
        mLambda(:, :, :, g1) = SSta.mLambdaFilt;
    end
    Time = toc;
    finaltime(k,1)=Time;
end
save('output\TimeExercisesTable3CPU.mat', 'finaltime');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% GPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
addpath('../DeCo_Package/codeGPU');
finaltime=NaN(2,1);
for kex=1:2
    
    bol.bias=kex;    
    % _______________________Load data and predictive ________
    load('input/UnBiasedPred.mat');
    vY=data.Y;
    if bol.bias==1
        mX=permute(data.fcst,[1 3 2]);
    else
        mX=permute(data.fcst(:,2:end,:),[1 3 2]);
    end
        % _____________________ Parameter setting __________________________
    iL=1;
    [cT, iM, iKL] = size(mX);
    iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
    vY = vY(end - cT + 1:end, iSelectSeries);
    %                       Parameter setting                                 %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    iDraws=100;
    Setting.cT = cT;                                    %%% Time dimension
    Setting.cN = 1000; %5000                            %%% Number of particles
    Setting.iDimOmega = iKL;                            %%% State space dimension
    Setting.dLambda = 0.95;                             %%% Lambda parameter
    Setting.iL = iL;                                    %%% L parameter
    Setting.iM = iM;                                    %%% M parameter
    Setting.iDraws = iDraws;                               %%% Draws from the predictive
    Setting.dKappa = 1;                                 %%% Resampling parameter
    Setting.iEstimate = 0.01;                           %%% Estimate or not the Lambda and Sigma
    Setting.dKappa = 1;                                 %%% Resampling parameter
    iLearning=0;
    iTau=9;
    if iLearning == 0                                   %%% Learning or not
        Setting.iTau = -1;
    else
        Setting.iTau = iTau;
    end
    %%-------------------- Setting the Sigma and Lambda matrix --------------%%
    mDataRand=[0.01,0.1];
    mSigmaTemp = [];
    vSigmaTemp = [];
    for i = 1:Setting.iL
        vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
        mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
    end
    Setting.Sigma = mSigmaTemp;
    Setting.vSigma = vSigmaTemp;
    % _____________________ Filtering __________________________________
    weightsCum = zeros(Setting.cT, Setting.iDimOmega, Setting.iDraws);
    %TrueweightsCumRange = zeros(Setting.cT, Setting.iDimOmega, 4);
    iCounter = 0;
    vVector = 1:Setting.iDraws:Setting.iM + 1;
    try
    tic;
    for i=1:size(vVector, 2) - 1
        disp(char('We are at the iteration: ', num2str(i)));
        mSigma = [];
        mXInit = zeros(Setting.cT * Setting.iDraws, Setting.iDimOmega * Setting.iDraws);
        mY = zeros(Setting.cT * Setting.iDraws, Setting.iDraws);
        iTemp = 0; iSeries = 0; iIter = 1;
        for g = (vVector(i):vVector(i + 1) - 1)
            %%% ---------------- Construct the mX and mError -------------- %%%
            mXInit(iTemp + 1:iTemp + Setting.cT, iSeries + 1:iSeries + Setting.iDimOmega) = reshape(mX(:, g, :), Setting.cT, Setting.iDimOmega);
            mY(iTemp + 1:iTemp + Setting.cT, iIter) = vY(:, 1);
            iTemp = iTemp + Setting.cT;
            iSeries = iSeries + Setting.iDimOmega;
            mSigma = blkdiag(mSigma, Setting.Sigma); %% This is the Sigma for the PF
            iIter = iIter + 1;
        end
        %%%%%%%% Load on the GPU %%%%%%%%
        mXGPU = gpuArray(mXInit);
        mYGPU = gpuArray(mY);
        %%%%%%%%%%% Initializzation of the elements %%%%%%%%%%%%%%%%%%%%%%%%%%%
        [S0_GPU] = initPSGPUSim(Setting, mYGPU, mXGPU);                    %% COMMENT: the initialization in BCRVD(2013) is done using values larger than the priors (and not equal to priors as in the empirical application and in the DeCo toolbox). To replicate BCRVD(2013) results we use the same values
        %%%%%%%%%%%%%%%% Clean the memory before Particle Filter %%%%%%%%%%%%%%
        clear mYGPU mXGPU;
        %%%%%%%%%%%%%%%%%%%%%%% Initialize the states %%%%%%%%%%%%%%%%%%%%%%%%%
        [SSta0_GPU] = initStaGPU(Setting);
        %%%%%%%%%%%%%%%% Preparing dataset for PF   %%%%%%%%%%%%%%%%%%%%%%%%%%%
        mXTot = zeros(Setting.iDimOmega * Setting.iDraws, Setting.cT, Setting.iL);
        mError = zeros(Setting.cT, Setting.iDimOmega * Setting.iDraws, Setting.iL);
        z = 1;
        for k = 1:Setting.iDraws
            mXTot(z:z + Setting.iDimOmega - 1, :) = reshape(mX(:, k + iCounter, :), Setting.cT, Setting.iDimOmega)';  
            vYTemp = kron(ones(1, Setting.iDimOmega), vY(:, 1));    %%%%Reshape the observations to create the learning errors
            mError(:, z:z + Setting.iDimOmega - 1) =  (vYTemp - reshape(mX(:, k + iCounter, :), Setting.cT, Setting.iDimOmega)).^2; %%%  Construct the errors
            z = z + Setting.iDimOmega;
        end
        %%%%%%%%%%%%%%%%%%%%% Starting the particle filter   %%%%%%%%%%%%%%%%%%
        SettingGPU = struct('cN', gpuArray(Setting.cN),'iDimOmega', gpuArray(Setting.iDimOmega),...
            'Sigma', gpuArray(Setting.Sigma),'vSigma',gpuArray(Setting.vSigma), ...
            'cT', gpuArray(Setting.cT), 'iDraws', gpuArray(Setting.iDraws), ...
            'iTau', gpuArray(Setting.iTau), 'dLambda', gpuArray(Setting.dLambda), ...
            'iL', gpuArray(Setting.iL), 'iM', gpuArray(Setting.iM), 'iEstimate', gpuArray(Setting.iEstimate), 'dKappa', gpuArray(Setting.dKappa));
        %%%%%%%%%%%%%%%%%%%%% Run the Particle filter on the GPU %%%%%%%%%%%%%%
        [~, SSta] = PFGPU(SettingGPU, S0_GPU, SSta0_GPU, gpuArray(vY), gpuArray(mXTot), gpuArray(mError));
        %%%%%%%%%%%%%%%%%% Go back to the CPU and do managing %%$$$$$$%%%%%%%%%
        mXFilt = gather(SSta.mXFilt);
        mYCum = gather(reshape(SSta.mYPred, Setting.cT * Setting.iDraws,  Setting.cN, []));
        clear SSta S S0_GPU SSta0_GPU;     %% Cancel old variables to freee memory in the GPU
        %%%%%%%%%% Managing the mean and quantile of the predictive %%%%%%%%%%%
        for g = 1:Setting.iDraws
            for kk=1:Setting.iL
                mYPred(:, 1, kk, g + iCounter) =  median(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk)), 2);
                mYPred(:, 2, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.025)';
                mYPred(:, 3, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.975)';
            end
        end
        iCounter = iCounter + Setting.iDraws;
        %%%%%%%%%%%%%%%%%%%% Managing the Weights %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for ga = 1:Setting.cT
            z = 1;
            for j = 1:Setting.iDraws
                weightsCum(ga, :, j + vVector(i) - 1) = logmul(squeeze(mXFilt(z:z + Setting.iDimOmega - 1, ga))');
                z = z + Setting.iDimOmega;
            end
        end
    end
    catch
      disp('-----------------------------------------------------------------');
      disp(char('The dimension of the particles', num2str(Setting.cN), 'and the dimension of the blocks', num2str(iDraws)));
       disp(char('are too big for the RAM available on the GPU card.')); 
       disp(char('Please try to decrease the dimension of the particles or the blocks or both.'));
       disp('-----------------------------------------------------------------');
     break;
     end
    Time = toc;
    finaltime(kex,1)=Time;
end
save('output/TimeExercisesTable3GPU.mat', 'finaltime');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of Table 3 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    Reproduce Table 4
%%%%%%%%%%%%%%%%%%%% CPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
addpath('../DeCo_Package/codeCPU');
load(strcat('input/Total_medium.mat'));
dire='output/Empirical/';
[cT, iM, iL, iKL] = size(mX);
iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
vY = vY(end - cT + 1:end, iSelectSeries);
% Parameter setting
Setting.cN = 1000; %5000                            %%% Number of particles
Setting.iDimOmega = iKL;                            %%% State space dimension
Setting.iM = iM;                                    %%% M parameter
Setting.iL = iL;                                    %%% L parameter
Setting.cT = cT;                                    %%% Time dimension
Setting.dLambda = 0.95;                             %%% Lambda matrix
Setting.iEstimate = 0.01;                           %%% Estimate or not the Lambda and Sigma
Setting.dKappa = 1;                                 %%% Resampling parameter
iLearning = 0;
iTau=9;
if iLearning == 0                                   %%% Learning or not
    Setting.iTau = -1;
else
    Setting.iTau = iTau;
end
%%-------------------- Setting the Sigma and Lambda matrix --------------%%
mDataRand = [0.1,0.25;0.01,0.4];
mSigmaTemp = [];
vSigmaTemp = [];
for i = 1:Setting.iL
   vSigmaTemp = [vSigmaTemp; log(mDataRand(i, 1))];
   mSigmaTemp =  [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega)]; 
end
Setting.Sigma = mSigmaTemp;
Setting.vSigma = vSigmaTemp;
%%% ------------ Create the matrix that will save the results ----------%%% 
mWeights = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
mYPred = zeros(Setting.cT, 3,  Setting.iL, Setting.iM);
mSigma = zeros(Setting.cT, Setting.iL, Setting.iM);
mLambda = zeros(Setting.cT, Setting.iDimOmega, Setting.iL, Setting.iM);
%%% ------------------ Start run filtering and timing ----------------- %%%
tic;
parfor g4 = 1:Setting.iM
    mError = zeros(Setting.cT, Setting.iDimOmega, Setting.iL);
    mXInit = zeros(Setting.cT, Setting.iDimOmega, Setting.iL);
    disp(char('We are at the iteration: ', num2str(g4)));
    %%% ---------------- Construct the mX and mError ------------------ %%%
    for i = 1:Setting.iL
        mXInit(:, :, i) = reshape(mX(:, g4, i, :), Setting.cT, Setting.iDimOmega);
        mError(:, :, i) = (Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, i)) - mXInit(:, :, i)).^2;
    end
    [S0] = initPSCPU(Setting, vY, mXInit);
    [SSta0] = initStaCPU(Setting);
    [S, SSta] = PFCPU(Setting, S0, SSta0, vY, mXInit, mError);
    mWeightsTemp = zeros(Setting.cT, Setting.iDimOmega, Setting.iL);
    for i = 1:Setting.cT
        for z = 1:Setting.iL
             mWeightsTemp(i, :, z) = logmul(SSta.mXFilt(i, :, z));  
        end     
    end    
    mWeights(:, :, :, g4) = mWeightsTemp; 
    mYPred(:, :, :, g4) = SSta.mYPred;
    mSigma(:, :, g4) = SSta.mSigmaFilt;
    mLambda(:, :, :, g4) = SSta.mLambdaFilt;
end
Time=toc;
save('output/TimeExercisesTable4CPU.mat', 'Time');
% Saving forecasts
save(strcat(dire, 'PredTotEmpCPU.mat'), 'mYPred');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% GPU  %%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
addpath('../DeCo_Package/codeGPU');
load('input/Total_medium.mat');
dire='output/Empirical/';
[cT, iM, iL, iKL] = size(mX);
iSelectSeries = 1:1:iL;             %%------ Select the iL series -------%%
vY = vY(end - cT + 1:end, iSelectSeries);
%%_____________________ Parameter setting __________________________
iDraws = 100;
Setting.cT = cT;                                    %%% Time dimension
Setting.cN = 1000; %5000                            %%% Number of particles
Setting.iDimOmega = iKL;                            %%% State space dimension
Setting.iM = iM;                                    %%% M parameter
Setting.iL = iL;                                    %%% L parameter
Setting.dLambda = 0.95;                             %%% Lambda matrix
Setting.iDraws = iDraws;
Setting.iEstimate = 0.01;                           %%% Estimate or not the Lambda and Sigma
Setting.dKappa = 1;                                 %%% Resampling parameter
iLearning = 0;
iTau=9;
if iLearning == 0                                   %%% Learning or not
    Setting.iTau = -1;
else
    Setting.iTau = iTau;
end
%%-------------------- Setting the Sigma and Lambda matrix --------------%%
mDataRand=[0.1,0.25;0.01,0.4];
mSigmaTemp = [];
vSigmaTemp = [];
for i = 1:Setting.iL
   vSigmaTemp = [vSigmaTemp, log(mDataRand(i, 1))];
   mSigmaTemp = [mSigmaTemp; log(mDataRand(i, 2)) * ones(1, Setting.iDimOmega);];
end
Setting.Sigma = mSigmaTemp;
Setting.vSigma = vSigmaTemp;
%%% ------------ Create the matrix that will save the results ----------%%% 
mWeightsCum = zeros(Setting.cT, Setting.iDimOmega, Setting.iDraws, Setting.iL);
mYPred = zeros(Setting.cT, 3, Setting.iL, Setting.iM);
TrueweightsCumRange = zeros(Setting.cT, Setting.iDimOmega, 4, Setting.iL);
vVector = 1:Setting.iDraws:Setting.iM + 1;
%%%%%%%-------- Run the particle filter for all the simulations ----- %%%%%
iCounter = 0;
try           %%%--------- Catch possible error on GPU memory ---------%%%
tic;
for i=1:size(vVector, 2) - 1
    disp(char('We are at the iteration: ', num2str(i))); 
    mSigma = []; 
    mXInit = zeros(Setting.cT * Setting.iDraws, Setting.iDimOmega * iDraws, Setting.iL);
    mY = zeros(Setting.cT * Setting.iDraws, Setting.iDraws, Setting.iL);
    iTemp = 0; iSeries = 0; iIter = 1;
    for g = (vVector(i):vVector(i + 1) - 1)
        %%% ---------------- Construct the mX and mError -------------- %%%
        for j = 1:iL
             mXInit(iTemp + 1:iTemp + Setting.cT, iSeries + 1:iSeries + Setting.iDimOmega, j) = reshape(mX(:, g, j, :), Setting.cT, Setting.iDimOmega); 
             mY(iTemp + 1:iTemp + Setting.cT, iIter, j) = vY(:, j);
        end
        iTemp = iTemp + Setting.cT;
        iSeries = iSeries + Setting.iDimOmega; 
        mSigma = blkdiag(mSigma, Setting.Sigma); %% This is the Sigma for the PF
        iIter = iIter + 1;
    end
    %%%%%%%% Load on the GPU %%%%%%%%
    mXGPU = gpuArray(mXInit);
    mYGPU = gpuArray(mY);
    %%%%%%%%%%% Initializzation of the elements %%%%%%%%%%%%%%%%%%%%%%%%%%%
    [S0_GPU] = initPSGPU(Setting, mYGPU, mXGPU);
    %%%%%%%%%%%%%%%% Clean the memory before Particle Filter %%%%%%%%%%%%%%
    clear mYGPU mXGPU;
    %%%%%%%%%%%%%%%%%%%%%%% Initialize the states %%%%%%%%%%%%%%%%%%%%%%%%%
    [SSta0_GPU] = initStaGPU(Setting);
    %%%%%%%%%%%%%%%% Preparing dataset for PF   %%%%%%%%%%%%%%%%%%%%%%%%%%%
    mXTot = zeros(Setting.iDimOmega * Setting.iDraws, Setting.cT, Setting.iL);
    mError = zeros(Setting.cT, Setting.iDimOmega * Setting.iDraws, Setting.iL);
    z = 1;
    for k = 1:Setting.iDraws
        for kk = 1:Setting.iL
            mXTot(z:z + Setting.iDimOmega - 1, :, kk) = reshape(mX(:, k + iCounter, kk, :), Setting.cT, Setting.iDimOmega)';   
            %vYTemp = kron(ones(1, Setting.iDimOmega), vY(:, kk));    %%%%Reshape the observations to create the learning errors
            vYTemp = Kronbsxfun(ones(1, Setting.iDimOmega), vY(:, kk));    %%%%Reshape the observations to create the learning errors
            mError(:, z:z + Setting.iDimOmega - 1, kk) =  (vYTemp - reshape(mX(:, k + iCounter, kk, :), Setting.cT, Setting.iDimOmega)).^2; %%%  Construct the errors 
        end
        z = z + Setting.iDimOmega;
    end
    %%%%%%%%%%%%%%%%%%%%% Starting the particle filter   %%%%%%%%%%%%%%%%%%
    SettingGPU = struct('cN', gpuArray(Setting.cN),'iDimOmega', gpuArray(Setting.iDimOmega),...
      'Sigma', gpuArray(Setting.Sigma),'vSigma',gpuArray(Setting.vSigma), ...
      'cT', gpuArray(Setting.cT), 'iDraws', gpuArray(Setting.iDraws), ...
      'iTau', gpuArray(Setting.iTau), 'dLambda', gpuArray(Setting.dLambda), ...
      'iL', gpuArray(Setting.iL), 'iM', gpuArray(Setting.iM), 'iEstimate', gpuArray(Setting.iEstimate), 'dKappa', gpuArray(Setting.dKappa));
    %%%%%%%%%%%%%%%%%%%%% Run the Particle filter on the GPU %%%%%%%%%%%%%%
    [~, SSta] = PFGPU(SettingGPU, S0_GPU, SSta0_GPU, gpuArray(vY), gpuArray(mXTot), gpuArray(mError));
    %%%%%%%%%%%%%%%%%% Go back to the CPU and do managing %%$$$$$$%%%%%%%%%
    mXFilt = gather(SSta.mXFilt);
    mLambda = gather(SSta.mLambdaFilt);
    mYCum = gather(reshape(SSta.mYPred, Setting.cT * Setting.iDraws,  Setting.cN, []));
    mSigma = gather(SSta.mSigmaFilt);
    clear SSta S S0_GPU SSta0_GPU;     %% Cancel old variables to freee memory in the GPU
    %%%%%%%%%% Managing the mean and quantile of the predictive %%%%%%%%%%% 
    for g = 1:Setting.iDraws
        for kk=1:Setting.iL
             mYPred(:, 1, kk, g + iCounter) =  median(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk)), 2);
             mYPred(:, 2, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.025)';
             mYPred(:, 3, kk, g + iCounter) =  quantile(squeeze(mYCum((g - 1) * Setting.cT + 1:(g * Setting.cT), :, kk))', 0.975)'; 
        end
    end
    iCounter = iCounter + Setting.iDraws; 
    %%%%%%%%%%%%%%%%%%%% Managing the Weights %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for ga = 1:Setting.cT    
        z = 1;
        for j = 1:Setting.iDraws
           for k = 1:Setting.iL
               mWeightsCum(ga, :, j + vVector(i) - 1, k) = logmul(squeeze(mXFilt(z:z + Setting.iDimOmega - 1, ga, k))');
           end      
           z = z + Setting.iDimOmega;
        end 
    end 
end
catch
   disp('-----------------------------------------------------------------');
   disp(char('The dimension of the particles', num2str(Setting.cN), 'and the dimension of the blocks', num2str(iDraws)));
   disp(char('are too big for the RAM available on the GPU card.')); 
   disp(char('Please try to decrease the dimension of the particles or the blocks or both.'));
   disp('-----------------------------------------------------------------');
 break;
end
Time = toc;
save('output/TimeExercisesTable4GPU.mat', 'Time');
% Saving forecasts
save(strcat(dire, 'PredtotEmpGPU.mat'), 'mYPred');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of Table 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;