function[F]=h3pdcsd_gen_det(N,T,b,p,q,y,x,w,db,dnb,tr,phi)
% SOS: this code uses 0.15 for trimming in lines 33 and 34 because p is
% large in the application
supF=zeros(1,b+1);
qs=(b+2)*q;
I=eye(T);
bs=b+1;
R=[eye(q) -eye(q)];
if bs>1
    for d=2:bs
        R=[R zeros((d-1)*q,q); zeros(q,(d-1)*q) eye(q) -eye(q)];
    end
end
if b>0
    if p>0
        [finalbreaksb]=pcdpcsd_gen_det(N,T,y,x,w,db,dnb,b,p,q,tr);
    else
        [~,finalbreaksb]=dpcsd_gen_det(N,T,y,w,[db dnb],b,q,tr);
    end
end
% Step 3: Get the SSRs for the extra break
if b==0
    kstar=[0 T];
else
    kstar=[0; finalbreaksb'; T];
end
for m=2:b+2
    Tstar=kstar(m)-kstar(m-1);
    wstar=w(kstar(m-1)+1:kstar(m),:,:);
    ystar=y(kstar(m-1)+1:kstar(m),:,:);
    dbstar=db(kstar(m-1)+1:kstar(m),:);
    dnbstar=dnb(kstar(m-1)+1:kstar(m),:);
    if p>0
        xstar=x(kstar(m-1)+1:kstar(m),:,:);
        [finalbreaks]=pcdpcsd_gen_det(N,Tstar,ystar,xstar,wstar,dbstar,dnbstar,1,p,q,tr);
    else
        [~,finalbreaks]=dpcsd_gen_det(N,Tstar,ystar,wstar,[dbstar dnbstar],1,q,tr);
    end
    finalbreaks=finalbreaks+kstar(m-1);
    if b==0
        k=sort([finalbreaks'; T]);
    else
        k=sort([finalbreaks'; finalbreaksb'; T]);
    end
    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*bs+1:q*bs+q,j)=[zeros(k(bs,1),q);w(k(bs,1)+1:T,:,j)];
        if bs>1
           for i=2:bs
                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
    ndb=size(db,2);
    ndbs=(b+2)*ndb;
    D=zeros(T,ndbs,1);
    D(:,1:ndb,1)=[db(1:k(1),:,1);zeros(T-k(1),ndb)];
    D(:,ndb*bs+1:ndb*bs+ndb,1)=[zeros(k(bs,1),ndb);db(k(bs,1)+1:T,:,1)];
    if bs>1
       for i=2:bs
            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];
    if p>0
        WTtilde=zeros(T,qs,N);
        xtilde=zeros(T,p,N);
        ytilde=zeros(T,1,N);
        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
        deltadenom=WW-WX*(XX\WX');
        deltanum=WY-WX*(XX\XY);
        delta=deltadenom\deltanum;
        if phi>0
            for i=1:N
                epshat(:,i)=ytilde(:,:,i)-xtilde(:,:,i)*(XX\XY)-WTtilde(:,:,i)*delta+xtilde(:,:,i)*(XX\WX')*delta;
                htilda(:,:,i)=WTtilde(:,:,i)-xtilde(:,:,i)*(XX\WX');
                for t=1:T
                    h(t,:,i)=htilda(t,:,i)*epshat(t,i);
                end
            end
            pt=floor(4*(T/100)^(2/9));
            Psij=zeros(qs,qs,pt+1);
            for j=1:pt+1
                for i=1:N
                    for t=(j-1)+1:T
                        Psij(:,:,j)=Psij(:,:,j)+h(t,:,i)'*h(t-(j-1),:,i)/(N*T);
                    end
                end
            end
            PsiV=Psij(:,:,1);
            for j=2:pt+1
                PsiV=PsiV+(1-(j-1)/(pt+1))*(Psij(:,:,j)+Psij(:,:,j)');
            end
            Vartheta=((deltadenom\PsiV)/deltadenom)*(N*T); 
            supF(1,m-1)=(delta'*R((m-2)*q+1:(m-1)*q,:)'*((R((m-2)*q+1:(m-1)*q,:)*Vartheta*R((m-2)*q+1:(m-1)*q,:)')\(R((m-2)*q+1:(m-1)*q,:)*delta)))/(q);
        else
            SSR=YY-XY'*(XX\XY)-(WY'-XY'*(XX\WX'))*(deltadenom\(WY'-XY'*(XX\WX'))');
            sigma2=SSR/(N*(T-p-qs)-p-qs);
            supF(1,m-1)=(delta'*R((m-2)*q+1:(m-1)*q,:)'*((R((m-2)*q+1:(m-1)*q,:)*(deltadenom\R((m-2)*q+1:(m-1)*q,:)'))\(R((m-2)*q+1:(m-1)*q,:)*delta)))/(q*sigma2);
        end
    else
        WTtilde=zeros(T,qs,N);
        ytilde=zeros(T,1,N);
        Zbar=[D mean(WT,3)];
        Q=I-Zbar*((Zbar'*Zbar)\Zbar');
        for i=1:N
            WTtilde(:,:,i)=Q*WT(:,:,i);
            ytilde(:,:,i)=Q*y(:,:,i);
        end
        WW=zeros(qs,qs);
        WY=zeros(qs,1);
        YY=0;
        for i=1:N
            WW=WW+WTtilde(:,:,i)'*WTtilde(:,:,i);
            WY=WY+WTtilde(:,:,i)'*ytilde(:,:,i);
            YY=YY+ytilde(:,:,i)'*ytilde(:,:,i);
        end
        deltadenom=WW;
        deltanum=WY;
        delta=deltadenom\deltanum;
        if phi>0
            for i=1:N
                epshat(:,i)=ytilde(:,:,i)-WTtilde(:,:,i)*delta;
                htilda(:,:,i)=WTtilde(:,:,i);
                for t=1:T
                    h(t,:,i)=htilda(t,:,i)*epshat(t,i);
                end
            end
            pt=floor(4*(T/100)^(2/9));
            Psij=zeros(qs,qs,pt+1);
            for j=1:pt+1
                for i=1:N
                    for t=(j-1)+1:T
                        Psij(:,:,j)=Psij(:,:,j)+h(t,:,i)'*h(t-(j-1),:,i)/(N*T);
                    end
                end
            end
            PsiV=Psij(:,:,1);
            for j=2:pt+1
                PsiV=PsiV+(1-(j-1)/(pt+1))*(Psij(:,:,j)+Psij(:,:,j)');
            end
            Vartheta=((deltadenom\PsiV)/deltadenom)*(N*T); 
            supF(1,m-1)=(delta'*R((m-2)*q+1:(m-1)*q,:)'*((R((m-2)*q+1:(m-1)*q,:)*Vartheta*R((m-2)*q+1:(m-1)*q,:)')\(R((m-2)*q+1:(m-1)*q,:)*delta)))/(q);
        else
            SSR=YY-(WY'*(deltadenom\WY));
            sigma2=SSR/(N*(T-qs)-qs);
            supF(1,m-1)=(delta'*R((m-2)*q+1:(m-1)*q,:)'*((R((m-2)*q+1:(m-1)*q,:)*(deltadenom\R((m-2)*q+1:(m-1)*q,:)'))\(R((m-2)*q+1:(m-1)*q,:)*delta)))/(q*sigma2);
        end
    end
end
[F,~]=max(supF);