/* SimpleDoubleSim1A.sas */
title 'Simple Double Measurement Regression Simulation';

data raw;
     /************ Set n and parameters ********************/
     n = 1000; mux = 4; phix = 2; beta0 = 1; beta1=0; psi = 1;
     omega11 = 1; omega22 = 3; corr = 0.5; /* corr(e1,e2) */
     /******************************************************/
     omega12 = corr*sqrt(omega11*omega22);
     /* The put statements write on the log file. */
     put ' *************************************************';
     put '            True omega12 = ' omega12 ;
     put ' *************************************************';
     /* Produce correlated errors by e1 = A1+B, e2 = A2+B
        with A1 and A2 independent, Var(A1)=va1, Var(A2)=va2,
        Cov(e1,e2) = E(A1+B)(A2+B)=Var(B)=omega12.
        Var(e1) = Var(A1+B) = va1+omega12 =>                */
     va1 = omega11-omega12; va2 = omega22-omega12;     
     do i = 1 to n;
        /* Latent model */
        X = mux + sqrt(phix)*rannor(378); 
        epsilon = sqrt(psi)*rannor(378);
        Y = beta0 + beta1*X + epsilon;
        /* Measurement model */
        A1 = sqrt(va1)*rannor(378); A2 = sqrt(va2)*rannor(378);
        B  = sqrt(omega12)*rannor(378); e1 = A1+B; e2 = A2+B; 
        W1 = X + e1; V1 = Y + e2;
        /* New A1 etc. to generate e3 and e4 */
        A1 = sqrt(va1)*rannor(378); A2 = sqrt(va2)*rannor(378);
        B  = sqrt(omega12)*rannor(378); e3 = A1+B; e4 = A2+B; 
        W2 = X + e3; V2 = Y + e4;
        output; /* Create a case */
     end;
     keep W1 W2 V1 V2;

proc calis pcorr pshort vardef=n nostand; 
 	 title2 'Fitted model incorrectly has omega12=0';
     var W1 W2 V1 V2;  /* Declare observed variables */
     lineqs        /* Model equations, separated by commas. */
        FY  = beta1*FX + epsilon,
        W1 = FX + e1, V1 = FY + e2,
        W2 = FX + e3, V2 = FY + e4;
     variance      /* Declare variance parameters. */
        FX = phix, epsilon = psi, 
        e1=omega11, e2=omega22, e3 = omega11, e4 = omega22;