%%This file is to generate data in a specific manner
% The inputs are the specifications for dimenesions/model/type
% The output is a 3-D cube, Variables-Time-Repetitions
function [Y, true] = DGP(dimension,Type,density)
%%
N = dimension(1);     % Number of Variables
T = dimension(2);     % Number of Timepoints
R = dimension(3);     % Number of Repetitions

Y = zeros(N,T,R);
A = zeros(N,N);

%% Generate coefficent matrix and covariance matrix for different DGP type
S = eye(N);
switch Type
    case 'Diagonal'
        A = 0.5*eye(N);
        Struct = eye(N);
        
    case 'IID'
        A = zeros(N);
        Struct = eye(N);
        
    case 'Diminishing-Diagonal'
        a = zeros(1,2*N-1);
        for n = 1:2*N-1
            a(n) = 0.3^(abs(N-n)+1);     %creates assistant vector for the VAR matrix
        end
        for n = 1:N
            A(n,:) = a(N-n+1:2*N-n);       %puts the vectors together for the VAR Matrix
        end
        Struct = eye(N);
        
    case 'Random'
        eigenModo = 2;
        while eigenModo>=0.99
            a_killer = rand(N,N)<(density);
            A = real(a_killer);
            A = (0.5)*A/max(abs(eig(A)));
            eigenModo = max(abs(eig(A)));
        end
        
        S_killer = rand(N,N)<(density);
        S_killer = triu(S_killer,1);
        B = real(S_killer)*density;
        Struct = (eye(N)-B)\diag(linspace(1,2,N)).^0.5;

        S =  Struct*Struct';
        
    case 'Network'
        noBlocks = 1/density;
        blockSize = N/noBlocks;
        blockSize = floor(blockSize);
        eigenModo = 2;
        while eigenModo>=0.99
            Sp = [];
            for n = 1:noBlocks
                a = eye(blockSize)*0.5+0.5*ones(blockSize,blockSize);
                Sp = [Sp zeros(size(Sp,1),size(a,2)); zeros(size(a,1),size(Sp,2)) a];%creates assistant vector for the VAR matrix
            end
            a = eye(N-size(Sp,1))*0.5+0.5*ones(N-size(Sp,1),N-size(Sp,2));
            Sp = [Sp zeros(size(Sp,1),size(a,2)); zeros(size(a,1),size(Sp,2)) a];%creates assistant vector for the VAR matrix
            A = Sp;
            A = (0.5)*A/max(abs(eig(A)));
            eigenModo = max(abs(eig(A)));
        end
        
        Struct = (eye(N)-A)\diag(linspace(1,2,N)).^0.5;

        S =  Struct*Struct';
end
%% Generate the output/observation data
for r = 1:R
    
    res = mvnrnd(zeros(N,1),S)';
    Ya(:,1) = res;
    for t = 2:T
        res = mvnrnd(zeros(N,1),S)';
        Ya(:,t) = A*Ya(:,t-1)+res;
    end
    Y(:,:,r) = Ya;
    
end

%% write results in a structure
true.A = full(A);
true.Sig = full(S);
true.Struct = full(Struct);
end
