function val = jsz_rotation(W, K1Q_X, K0Q_X, dt, BX, AX)
%     ## Inputs:
%     ##   W          : N*J,      vector of portfolio weights to fit without error.
%     ##   K1Q.X      : N*N
%     ##   K0Q.X      : N*1
%     ##   dt         : scalar,   length of period in years
%     ##   BX         : N*J  (BX, AX) are optional (saves time)
%     ##   AX         : 1*J
%     ##
%     ## Returns:  [K0Q.cP, K1Q.cP, rho0.cP, rho1.cP]
%     ##   K0Q.cP : N*1
%     ##   K1Q.cP : N*N
%     ##   rho0.cP : scalar
%     ##   rho1.cP : N*1
%     ##
%     ## r(t) = rho0.cP + rho1.cP'*cPt
%     ##      = 1'*Xt
%     ##      = 1 period discount rate (annualized)
%     ##
%     ## Under Q:
%     ##   X(t+1) - X(t)   = K0Q.X  + K1Q.X*X(t)  + eps_X(t+1),   cov(eps_X(t+1)) = Sigma_X
%     ##   cP(t+1) - cP(t) = K0Q.cP + K1Q.cP*X(t) + eps_cP(t+1),  cov(eps_cP(t+1)) = Sigma.cP
%     ## Where Sigma_X is chosen to match Sigma.cP
%     ##
%     ## cPt = W*yt  (cPt N*1, W is N*J)
%     ##     = W*AX' + W*BX'*Xt
%     ##     = WAXp + WBXp*Xt
%     ##
%     ## Delta(cP) = WBXp*Delta(Xt)
%     ##           = WBXp*(K1Q.X*Xt + sqrt(Sigma_X)*eps(t+1))
%     ##           = WBXp*(K1Q.X)*(WBXp\(cPt - WAXp)) + sqrt(Sigma.cP)*eps(t+1)
%     ##           = WBXp*(K1Q.X)/WBXp*cPt - WBXp*(K1Q.X)/WBXp*WAXp] + sqrt(Sigma.cP)*eps(t+1)
%     ##
%     ## rt = 1'*Xt  [annualized 1-period rate]
%     ##    = 1'*(WBXp\(cPt - WAXp))
%     ##    = [- 1'*(WBXp\WAXp)] + ((WBXp)'1)'*cPt
    
N = size(K1Q_X,1) ;
WBXp = W * BX' ;
WAXp = W * AX' ;
WBXpinv = inv(WBXp) ;

K1Q_cP = WBXp * K1Q_X * WBXpinv ;
K0Q_cP = WBXp * K0Q_X - K1Q_cP * WAXp ;

rho1_cP = WBXpinv' * ones(N,1) ;
rho0_cP = - WAXp' * rho1_cP    ;

val.K0Q_cP = K0Q_cP ;
val.K1Q_cP = K1Q_cP ;
val.rho0_cP = rho0_cP ;
val.rho1_cP = rho1_cP ;


end

 
    
    
    
    