function[finalbreaks]=pcdpcsd_gen_det(N,T,y,x,w,db,dnb,b,p,q,tr)
% Step 1: Get the initial beta
qs=(b+1)*q;
c=zeros(T,p+q,N);
for s=1:N
    c(:,:,s)=[x(:,:,s) w(:,:,s)];
end
[~,truebreaks]=dpcsd_gen_det(N,T,y,c,[db dnb],b,p+q,tr);
k=[truebreaks'; T];
% Step 2: regression with known breaks in c
CT=zeros(T,(b+1)*(p+q),N);
for j=1:N
    CT(:,1:p+q,j)=[[x(1:k(1),:,j) w(1:k(1),:,j)];zeros(T-k(1),p+q)];
    CT(:,(q+p)*b+1:(q+p)*b+(q+p),j)=[zeros(k(b,1),p+q);[x(k(b,1)+1:T,:,j) w(k(b,1)+1:T,:,j)]];
    if b>1
       for i=2:b
            CT(:,(q+p)*(i-1)+1:(q+p)*i,j)=[zeros(k(i-1,1),p+q);[x(k(i-1,1)+1:k(i,1),:,j) w(k(i-1,1)+1:k(i,1),:,j)];zeros(T-k(i,1),p+q)];
        end
    end
end
ndb=size(db,2)+size(dnb,2);
dbT=[db dnb];
ndbs=(b+1)*ndb;
D=zeros(T,ndbs,1);
D(:,1:ndb,1)=[dbT(1:k(1),:,1);zeros(T-k(1),ndb)];
D(:,ndb*b+1:ndb*b+ndb,1)=[zeros(k(b,1),ndb);dbT(k(b,1)+1:T,:,1)];
if b>1
   for i=2:b
        D(:,ndb*(i-1)+1:ndb*i,1)=[zeros(k(i-1,1),ndb);dbT(k(i-1,1)+1:k(i,1),:,1);zeros(T-k(i,1),ndb)];
    end
end
CTtilde=zeros(T,(b+1)*(p+q),N);
ytilde=zeros(T,1,N);
I=eye(T);
Zbar=[D mean(CT,3)];
Q=I-Zbar*((Zbar'*Zbar)\Zbar');
for s=1:N
    CTtilde(:,:,s)=Q*CT(:,:,s);
    ytilde(:,:,s)=Q*y(:,:,s);
end
CC=zeros((b+1)*(p+q),(b+1)*(p+q));
CY=zeros((b+1)*(p+q),1);
for s=1:N
    CC=CC+CTtilde(:,:,s)'*CTtilde(:,:,s);
    CY=CY+CTtilde(:,:,s)'*ytilde(:,:,s);
end
theta=CC\CY;
% Step 3: Now I need a new dependent variable
WT=zeros(T,qs,N);
for j=1:N
    WT(:,1:q,j)=[w(1:k(1),:,j);zeros(T-k(1),q)];
    WT(:,q*b+1:q*b+q,j)=[zeros(k(b,1),q);w(k(b,1)+1:T,:,j)];
    if b>1
       for i=2:b
            WT(:,q*(i-1)+1:q*i,j)=[zeros(k(i-1,1),q);w(k(i-1,1)+1:k(i,1),:,j);zeros(T-k(i,1),q)];
        end
    end
end
splitcoef=reshape(theta,p+q,b+1);
delta=splitcoef(p+1:p+q,:);
YT=zeros(T,1,N);
for i=1:N
    YT(:,1,i)=y(:,1,i)-WT(:,:,i)*delta(:);
end
xtilde=zeros(T,p,N);
YTtilde=zeros(T,1,N);
ndb=size(db,2);
ndbs=(b+1)*ndb;
D=zeros(T,ndbs,1);
D(:,1:ndb,1)=[db(1:k(1),:,1);zeros(T-k(1),ndb)];
D(:,ndb*b+1:ndb*b+ndb,1)=[zeros(k(b,1),ndb);db(k(b,1)+1:T,:,1)];
if b>1
   for i=2:b
        D(:,ndb*(i-1)+1:ndb*i,1)=[zeros(k(i-1,1),ndb);db(k(i-1,1)+1:k(i,1),:,1);zeros(T-k(i,1),ndb)];
    end
end
D=[D dnb];
Zbar=[D mean(x,3)];
Q=I-Zbar*((Zbar'*Zbar)\Zbar');
for s=1:N
    xtilde(:,:,s)=Q*x(:,:,s);
    YTtilde(:,:,s)=Q*YT(:,:,s);
end
XX=zeros(p,p);
XY=zeros(p,1);
for s=1:N
    XX=XX+xtilde(:,:,s)'*xtilde(:,:,s);
    XY=XY+xtilde(:,:,s)'*YTtilde(:,:,s);
end
beta0=XX\XY;
% Step 4:
SSRf=zeros(1,10);
beta=zeros(p,10);
count=0;
error=1;
while error>0
    count=count+1;
    % Step 4.1: keeping beta fixed, estimate delta and break dates
    YN=zeros(T,1,N);
    for i=1:N
        YN(:,1,i)=y(:,1,i)-x(:,:,i)*beta0(:,count);
    end
    [~,truebreaks]=dpcsd_gen_det(N,T,YN,w,[db dnb],b,q,tr);
    k=[truebreaks'; T];
    % Step 5: with the new breaks, estimate beta and delta and get the SSR
    WT=zeros(T,qs,N);
    for j=1:N
        WT(:,1:q,j)=[w(1:k(1),:,j);zeros(T-k(1),q)];
        WT(:,q*b+1:q*b+q,j)=[zeros(k(b,1),q);w(k(b,1)+1:T,:,j)];
        if b>1
           for i=2:b
                WT(:,q*(i-1)+1:q*i,j)=[zeros(k(i-1,1),q);w(k(i-1,1)+1:k(i,1),:,j);zeros(T-k(i,1),q)];
            end
        end
    end
    WTtilde=zeros(T,qs,N);
    xtilde=zeros(T,p,N);
    ytilde=zeros(T,1,N);
    D=zeros(T,ndbs,1);
    D(:,1:ndb,1)=[db(1:k(1),:,1);zeros(T-k(1),ndb)];
    D(:,ndb*b+1:ndb*b+ndb,1)=[zeros(k(b,1),ndb);db(k(b,1)+1:T,:,1)];
    if b>1
       for i=2:b
            D(:,ndb*(i-1)+1:ndb*i,1)=[zeros(k(i-1,1),ndb);db(k(i-1,1)+1:k(i,1),:,1);zeros(T-k(i,1),ndb)];
        end
    end
    D=[D dnb];
    Zbar=[D mean(x,3) mean(WT,3)];
    Q=I-Zbar*((Zbar'*Zbar)\Zbar');
    for i=1:N
        xtilde(:,:,i)=Q*x(:,:,i);
        WTtilde(:,:,i)=Q*WT(:,:,i);
        ytilde(:,:,i)=Q*y(:,:,i);
    end
    WW=zeros(qs,qs);
    WX=zeros(qs,p);
    XX=zeros(p,p);
    WY=zeros(qs,1);
    XY=zeros(p,1);
    YY=0;
    for i=1:N
        WW=WW+WTtilde(:,:,i)'*WTtilde(:,:,i);
        WX=WX+WTtilde(:,:,i)'*xtilde(:,:,i);
        XX=XX+xtilde(:,:,i)'*xtilde(:,:,i);
        WY=WY+WTtilde(:,:,i)'*ytilde(:,:,i);
        XY=XY+xtilde(:,:,i)'*ytilde(:,:,i);
        YY=YY+ytilde(:,:,i)'*ytilde(:,:,i);
    end
    betadenom=XX-WX'*(WW\WX);
    betanum=XY-WX'*(WW\WY);
    beta0(:,count+1)=betadenom\betanum;
    SSRf(1,count+1)=YY-WY'*(WW\WY)-(XY-WX'*(WW\WY))'*(betadenom\(XY-WX'*(WW\WY)));
    error=SSRf(1,count+1)-SSRf(1,count);
end
finalbreaks=k(1:b,1)';
% count