ols.src

来自「没有说明」· SRC 代码 · 共 689 行 · 第 1/2 页

        if __miss == 0 and ismiss(dta);
            errorlog "missing data found - __MISS set to 1 (listwise)";
            __miss = 1;
        endif;
        if __miss == 1;
            dta = packr(dta);
            if scalmiss(dta);
                goto errout(31);
            endif;
            nobs = rows(dta);
            mobs = tobs-nobs;
        endif;
        mn = meanc(dta);
        m = moment(dta,0)/nobs;
    else;
        clear mn,m,nc;
        do until eof(fin);
            y0 = readr(fin,nr);
            y0 = y0[.,vardx];
            if __miss == 1;
                y0 = packr(y0);
                nc = nc+rows(y0);
            elseif ismiss(y0);
                errorlog "missing data found - __MISS set to 1 (listwise)";
                __miss = 1;
                y0 = packr(y0);
                nc = nc+rows(y0);
            endif;
            if not scalmiss(y0);
                m = m+moment(y0,0);
                mn = mn + sumc(y0);
            endif;
        endo;
        if __miss == 1;
            if nc == 0;
                goto errout(31);
            endif;
            nobs = nc;
            mobs = tobs-nobs;
        endif;
        mn = mn/nobs;
        m = m/nobs;
    endif;

    if __miss == 2;
        constflg = indexcat(dotfeq(diag(m),diag(mn)^2),1);
    else;
        constflg = indexcat(dotfeq(diag(m),mn^2),1);
    endif;

    if scalmiss(constflg);
        constflg = 0;
    elseif rows(constflg) > 1;
        goto errout(35);
    endif;

    if constflg;
        cvec = packr(miss(seqa(1,1,rows(mn)),constflg));
        if __miss == 2;
            constvlu = mn[constflg,constflg];
            mn = mn[cvec,cvec];
        else;
            constvlu = mn[constflg];
            mn = mn[cvec];
        endif;
        m = m[cvec,cvec];
        nvar1 = rows(cvec);
        nvar = nvar1 - 1;
        if not idat;
            cnstname = indvars[constflg];
            indvars = indvars[packr(miss(seqa(1,1,rows(indvars)),constflg))];
        else;
            indvars = 0$+"X"$+ftocv(seqa(1,1,nvar),
                __vpad*(floor(log(nvar))+1),0);
        endif;
    endif;

    if __miss == 2;
        mn = diag(mn);
    endif;

    if const == 1 and constflg;
        const = 0;
    endif;

    if const or constflg;
      cov = m - mn*mn';
    else;
      cov = m;
    endif;

    k = diag(cov);
    cyy = k[nvar1];
    std = sqrt(k);
    cxy = cov[1:nvar,nvar1];
    cxx = cov[1:nvar,1:nvar];
    cor = cov./std./std';

    oldtrp = trapchk(1);
    trap 1,1;
    cxxi = invpd(cxx);
    trap oldtrp,1;
    if scalmiss(cxxi);
        goto errout(30);
    endif;

    b = cxxi*cxy;

    if const or constflg;
        constant = (mn[nvar1]-mn[1:nvar]'*b)/constvlu;
        if constflg == 0;
            indvars = "CONSTANT"|indvars;
        endif;
    endif;
    if idat and __altnam $/= 0;
        vnames = __altnam;
        if const == 0 and constflg;
            k = packr(miss(seqa(1,1,rows(vnames)),constflg));
            vnames = vnames[constflg]|vnames[k];
        endif;
        depvar = vnames[rows(vnames)];
        indvars = vnames[1:rows(vnames)-1];
    else;
        vnames = indvars|depvar;
    endif;

    if rows(indvars) == nvar and (const or constflg);
        if not idat;
            indvars = cnstname|indvars;
        else;
            indvars = "CONSTANT"|indvars;
        endif;
    endif;

    if const or constflg;
       df = nobs-nvar-1;
    else;
       df = nobs-nvar;
    endif;

    if df == 0;
        goto errout(32);
    elseif df<0;
        goto errout(31);
    endif;

    sse = cyy-b'*cxy;
    if const or constflg;
        k = -cxxi*mn[1:nvar]/constvlu;
        vc = (sse/df)*(((1/constvlu-mn[1:nvar]'*k)/constvlu|k)~(k'|cxxi));
        stderr = sqrt(diag(vc));
        t = (constant|b)./stderr;
        tv = nobs*cyy;
    else;
        vc = (sse/df)*cxxi;
        stderr = sqrt(diag(vc));
        t = b./stderr;
        tv = nobs*(cyy - mn[nvar1]^2);
    endif;
    sse = nobs*sse;
    rsq = (tv - sse)/tv;
    rbsq = 1-(1-rsq)*((nobs-1)/df);
    fstat = (rsq/(1-rsq))*(df/nvar);
    if fstat>0;
        pvf = cdffc(fstat,nvar,df);
    else;
        pvf = mss;
    endif;
    pvt = 2*cdftc(abs(t),df);

    if sse > 0;
      stdest = sqrt(sse/df);
    else;
      stdest = error(0);
    endif;

    stdb = b.*(std[1:nvar]/std[nvar1]);     /* Standardized coefficients  */
    if const or constflg;
        stdb = mss|stdb;
    endif;

    if _olsres;
        old = ndpcntrl(0,0);
        call ndpcntrl(1,1);
        if idat;
            if constflg;
                dta = dta[.,cvec];
            endif;
            u = dta[.,nvar1]-dta[.,1:nvar]*b - constant*constvlu;
            ndpclex;
            if __miss;
                u0 = packr(u);
            else;
                u0 = u;
            endif;
            dwstat = sumc((trimr(u0,1,0)-trimr(u0,0,1))^2)/(u0'*u0);
            clear u0;
        else;
            prcn = 8;
            if _olsres == 4;
                prcn = 4;
            endif;
            create fout = ^_olsrnam with u,1,prcn;
            if fout == -1;
                errorlog "Can't open temporary file for residuals";
                end;
            endif;
            call seekr(fin,1);
            clear dwstat,u2,i;
            do until eof(fin);
                i = i + 1;
                y0 = readr(fin,nr);
                y0 = y0[.,vardx];
                if __miss == 2;
                    y0 = missrv(y0,0);
                endif;
                if constflg;
                    y0 = y0[.,cvec];
                endif;
                u = y0[.,nvar1]-y0[.,1:nvar]*b - constant*constvlu;
                ndpclex;
                if writer(fout,u) /= rows(u);
                    errorlog "ERROR - disk full, Durbin-Watson statistic no"\
                        "t computed";
                    end;
                endif;
                if __miss;
                    u = packr(u);
                endif;
                u2 = u2+u'*u;
                if nr > 1;
                  dwstat = dwstat+sumc((trimr(u,1,0)-trimr(u,0,1))^2);
                endif;
                if i > 1 and i < tobs;
                  dwstat = dwstat + (u[1] - dd)^2;
                endif;
                dd = u[rows(u)];
            endo;
            dwstat = dwstat/u2;
        endif;
        call ndpcntrl(old,0xffff);
    else;
        u = 0;
        dwstat = 0;
    endif;
    if const or constflg;
        b = constant|b;
    endif;

    if __output;
        print ftos(nobs,"Valid cases:  %*.*lf",20,0);;
        print ftos(depvar,"      Dependent variable:%*.*s",20,8);

        print ftos(mobs,"Missing cases:%*.*lf",20,0);;
        print "      Deletion method:               ";;
        if __miss == 0;
            print "    None";
        elseif __miss == 2;
            print "Pairwise";
        else;
            print "Listwise";
        endif;

        print ftos(tv,"Total SS:     %*.*lf",20,3);;

        print ftos(df,"      Degrees of freedom:%*.*lf",20,0);

        print ftos(rsq,"R-squared:    %*.*lf",20,3);;
        print ftos(rbsq,"      Rbar-squared:      %*.*lf",20,3);
        print ftos(sse,"Residual SS:  %*.*lf",20,3);;
        print ftos(stdest,"      Std error of est:  %*.*lf",20,3);
        str = ftos(nvar,"F(%*.*lf,",1,0) $+ ftos(df,"%*.*lf):             "
            ,1,0);
        str = strsect(str,1,15) $+ ftos(fstat,"%*.*lf",19,3);
        print str;;
        print ftos(pvf,"      Probability of F:  %*.*lf",20,3);

        if _olsres;
            print ftos(dwstat,"Durbin-Watson:%*.*lf",20,3);
        endif;
        print;
        print "                         Standard                 Prob   Sta"\
            "ndardized  Cor with";
        print "Variable     Estimate      Error      t-value     >|t|     E"\
            "stimate    Dep Var";

        print "------------------------------------------------------------"\
            "-------------------";
        omat = indvars~b~stderr~t~pvt~stdb;
        if const or constflg;
            omat = omat~(mss|cor[1:nvar,nvar1]);
        else;
            omat = omat~cor[1:nvar,nvar1];
        endif;
        mask = 0~1~1~1~1~1~1;
        let fmt[7,3] = "-*.*s" 9 8 "*.*lf" 12 6 "*.*lf" 12 6 "*.*lf" 12 6 ""\
            "*.*lf" 10 3 "*.*lf" 12 6 "*.*lf" 12 6;
        ms = ftos(mss,"%*.*lf",1,0);
        msym "---  ";
        call printfm(omat,mask,fmt);
        msym ^ms;
    endif;
    if fin > 0;
        fin = close(fin);
    endif;
    if fout > 0;
        fout = close(fout);
    endif;

    m = (1~mn')|(mn~m);

    retp(vnames,nobs*m,b,stdb,vc,stderr,stdest,cor,rsq,u,dwstat);

ERROUT:
        pop be;
    if be == 34;
        errorlog "ERROR: File not found: " $+ dataset;
    elseif be == 30;
        errorlog "ERROR: covariance matrix of independent variables is sing"\
            "ular.";
    elseif be == 31;
        errorlog "ERROR: system underdetermined";
    elseif be == 32;
        errorlog "ERROR: same number columns as rows";
    elseif be == 33;
        errorlog "ERROR: too many missings";
    elseif be == 35;
        errorlog "ERROR: no variation in at least one independent variable";
    else;
        errorlog "Coefficients vector is an error code: " $+ ftos(be,"%*.*l"\
            "f",1,0);
    endif;
    if fin > 0;
        fin = close(fin);
    endif;
    if fout > 0;
        fout = close(fout);
    endif;
    retp(0,0,error(be),0,0,0,0,0,0,0,0);

endp;