function [yt,Xt,WN,So,Sm] = gendata(N,T,beta,s2m,s2n,rho1,rho2,type,cut,ss);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spatial weighting matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Filling the grid -  N observations
SQN=int32(N^0.5);
units=int8(zeros(2,N));
for i=1:SQN
    for j=1:SQN
        units(1,j+SQN*(i-1))=i;
        units(2,j+SQN*(i-1))=j;
    end
end

%spatial weighting matrix: Queen or distance
WN=sparse((zeros(1,1)));
if type==1 %rook
    for i=1:N
        for j=1:N
            if i==j
                WN(i,i)=0;
            else
                temp1=double(abs(units(1,i)-units(1,j)));
                temp2=double(abs(units(2,i)-units(2,j)));
                WN(i,j)=((temp1 == 1 & temp2 == 0) |  (temp1 == 0 & temp2 == 1));
            end
        end
    end
end
if type==2 % distance
    for i=1:N
        for j=1:N
            if i==j
                WN(i,i)=0;
            else
                temp1=double(abs(units(1,i)-units(1,j)));
                temp2=double(abs(units(2,i)-units(2,j)));
                dist=sqrt(temp1^2+ temp2^2);
                if dist < cut
                    WN(i,j)=1/dist;
                end
            end
        end
    end
end

clear temp1;
clear temp2;
clear units;
clear dist;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum row normalization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mmax=max(sum(WN,2)');
WN=sparse((inv(diag(mmax))*WN));
clear mmax

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate Zm, Xt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',15);
Zm=sparse(kron(ones(T,1),speye(N,N)));
xt= 5*Zm*(rand(N,1)-0.5)+2.5*(rand(N*T,1)-0.5);
Xt=[xt,ones(N*T,1)];
IN=speye(N);
clear xt Zm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%true model yT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%whos;
yt=Xt(:,1)*beta(1,1)+beta(2,1)*Xt(:,2);
mm=((IN-rho1*WN)\IN)*sqrt(s2m)*randn(N,1);

for t=1:T
    yt(((t-1)*N+1):((t-1)*N+N),1)=yt(((t-1)*N+1):((t-1)*N+N),1)+mm;
end
clear mm

for t=1:T
    yt(((t-1)*N+1):((t-1)*N+N),1)=yt(((t-1)*N+1):((t-1)*N+N),1)+...
        sparse((IN-rho2*WN)\IN)*(sqrt(s2n)*randn(N,1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Selector matrices to be provided 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',10);
skip=int8(zeros(N,1));

So=speye(2*N,2*N);
Sm=speye(2*N,2*N);

pack

for i=3:T
    Sto=speye(N,N);
    Stm=speye(N,N);
    nskip=round(ss*N);
    k=1;
    l=1;
    clear skip
    skip=int32(zeros(1,1));
    while size(skip,1) < nskip
        temp=int32(fix(rand(1,1)*N));
        if temp*int32(ones(size(skip,1),1))-skip & temp > 0
            skip(k,1)=int32(temp);
            k=k+1;
        end
    end

    sskip=sort(skip);
    smis=int32(zeros(1,1));
    smis=[];
    for i=1:N
        if i~=skip
            smis=[smis;i];
        end
    end
    ssmis=sort(smis);

    for j=1:nskip
        Sto(sskip(j,1),:)=[];
        sskip=sskip-1;
    end

    for j=1:(N-nskip)
        Stm(ssmis(j,1),:)=[];
        ssmis=ssmis-1;
    end
    
    ssio1=size(So,1);
    ssio2=size(So,2);
    ssim1=size(Sm,1);
    ssim2=size(Sm,2);

    sio1=size(Sto,1);
    sio2=size(Sto,2);
    sim1=size(Stm,1);
    sim2=size(Stm,2);
    
%     So = sparse(blkdiag(full(So),full(Sto)));
%     Sm = sparse(blkdiag(full(Sm),full(Stm)));

    %18 percent of every wave missing
    So=[So,zeros(ssio1,sio2)];
    So=[So;zeros(sio1,ssio2+sio2)];
    Sm=[Sm,zeros(ssim1,sim2)];
    Sm=[Sm;zeros(sim1,ssim2+sim2)];

    So(ssio1+1:ssio1+sio1,ssio2+1:ssio2+sio2)=Sto;
    Sm(ssim1+1:ssim1+sim1,ssim2+1:ssio2+sim2)=Stm;

end
Sm(1:2*N,:)=[];
clear smis ssmis skip sskip Stm Sto temp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
