%% this file is for the results similar to Lewbel 1991

clear all
close all
clc

H = 100;
G = 14;
T = 1;

START = 1997
END = 2006

count = 1;

mem = input('type of household:\n (1) 1 member\n (2) 2 members\n (3) 2-3 members\n (4) 2-4 members\n');
warning off;

for S = START:T:2006
    start_year = S;
    
    yg = zeros(H,G,T);
    te = zeros(H,T);
    w = [];
    
    for t = 1:T
        if mem == 1
            filename = ['data/' num2str(START) '_' num2str(END) '/dat' num2str(start_year+t-2) '_' num2str(start_year+t-1) 'def_mem1.csv']
        elseif mem == 2
            filename = ['data/' num2str(START) '_' num2str(END) '/dat' num2str(start_year+t-2) '_' num2str(start_year+t-1) 'def_mem2.csv']
        elseif mem == 3
            filename = ['data/' num2str(START) '_' num2str(END) '/dat' num2str(start_year+t-2) '_' num2str(start_year+t-1) 'def_mem2_3.csv']
        elseif mem == 4
            filename = ['data/' num2str(START) '_' num2str(END) '/dat' num2str(start_year+t-2) '_' num2str(start_year+t-1) 'def_mem2_4.csv']
        end
      
        data=csvread(filename,1,1);
        
        data(:,G+2) = [];
        
        yg(:,:,t) = data(:,2:G+1);
        te(:,t) = data(:,1);
        
        for g=1:G
            bs(:,g,t) =  yg(:,g,t);
        end
        
        clear data
    end
    
    w = reshape(bs,H,G*T);
    
    %% remove budget shares >1 and <0
    [row col]= find(w>=1 | w<0);
    if length(row)>0
        for i =1:length(row)
            w(row(i),col(i))= mean(w(:,col(i)));
        end
    end
    
    %% remove miscellaneous
    w(:,[1:T]*G)=[];
    
    %% center data
    [H N] = size(w);
    ss = std(w);
    mm = mean(w);
    ss1 = ones(size(w ,1) , 1)*ss;
    mm1 = ones(size(w, 1) , 1)*mm;
    x = w-mm1;
    
    %% Build functional forms
    FF(1,:) = te.*log(te);
    i=1;
    for a= [-3 -2 -1 1 2 3]
        Q1(i,:) = te.^a;
        Q2(i,:) = (log(te)).^a;
        i= i+1;
    end
    FF = [FF;Q1;Q2];
    FF = FF';
    clear Q1 Q2
    
    %% LDU
    G =13;
    for r=1:G-1
        C = (x(1:H,:))'*(FF(1:H,:))/(H-1);
        
        [evec eval] = eig(C);
        
        AUTO(:,S-1997+1) = diag(eval); 
        
        A = sparse(C);
        [L1,U1,P,Q] = lu(A);
        clear P Q
        L1 = full(L1);
        U1 = full(U1);
        D = diag(diag(U1));
        U1 = inv(D)*U1;
        [D1,IX] = sort((diag(D)));                % sort eigenvalues and eigenvectors
        D1 = flipud(D1);
        IX = flipud(IX);
        Q = eye(size(C,1));
        Q = Q(:,IX);
        
        LLL = L1*diag(D1)*U1 - L1*inv(Q')*inv(L1)*C*inv(U1)*inv(Q)*U1;
        
        P = L1*inv(Q')*inv(L1);
        Q = inv(U1)*inv(Q)*U1;
        L = L1;
        U = U1;
        D = diag(D1);
        
        clear L1 U1 D1
        
        L11 = L(1:r,1:r);
        L21 = L(r+1:end,1:r);
        L22 = L(r+1:end,r+1:end);
        
        D1 = D(1:r,1:r);
        D2 = D(r+1:end,r+1:end);
        
        U11 = U(1:r,1:r);
        U12 = U(1:r,r+1:end);
        U22 = U(r+1:end,r+1:end);
        
        dd2 = diag(D2);
        
        vD2 = vech(D2);
        
        Delta = zeros((G-r)^2,G-r);
        
        for j = 1:(G-r)
            Delta((j-1)*(G-r)+j,j) = 1;
        end
        
        Y = [-inv(L22)*L21*inv(L11) inv(L22)];
        K = [-inv(U11)*U12*inv(U22);inv(U22)];
        
        DUP = duplication(G);
        DUPD = inv(DUP'*DUP)*DUP';
        
        Vvect = 2*DUPD*kron(C,C)*DUPD'*H;
        
        V = 2*kron(C,C)*H;        
        
        
        W1 = Delta'*kron(K',Y)*kron(Q',P)*V*kron(Q,P')*kron(K,Y')*Delta;
        W2 = Delta'*kron(K',Y)*V*kron(K,Y')*Delta;

        STAT1(S-1997+1,r) = H*dd2'*inv(W1)*dd2;
        STAT2(S-1997+1,r) = H*dd2'*inv(W2)*dd2;

        dof(S-1997+1,r) = G-r;
        PV1(S-1997+1,r) = 1-chi2cdf(STAT1(S-1997+1,r),dof(S-1997+1,r));
        PV2(S-1997+1,r) = 1-chi2cdf(STAT2(S-1997+1,r),dof(S-1997+1,r));
    end
    G = 14;
    clear FF
end

for i = 1:size(PV1,1)
    n1 = find(PV1(i,:)>0.1);
    nfac190(i,1) = n1(1,1);
    n2 = find(PV1(i,:)>0.05);
    nfac195(i,1) = n2(1,1);
end

clear n1 n2

for i = 1:size(PV2,1)
    n1 = find(PV2(i,:)>0.1);
    nfac290(i,1) = n1(1,1);
    n2 = find(PV2(i,:)>0.05);
    nfac295(i,1) = n2(1,1);
end