Longitudinal Data: Variable Bandwidth Smoothing

/*--------------------------------------------------------------
                    SAS Sample Library

       Name: ssmex09.sas
       Description: Example program from SAS/ETS User's Guide,
              The SSM Procedure
       Title: Longitudinal Data: Variable Bandwidth Smoothing
       Product: SAS/ETS Software
        Keys: Longitudinal Data, Variable Bandwidth Smoothing
        PROC: SSM
        Notes:

--------------------------------------------------------------*/

 title;

 data Difficult;
 input x y;
 datalines;
 0.002 0.040
 0.011 0.009
 0.013 0.719
 0.016 0.199
 0.017 -0.409
 0.019 0.410
 0.021 -0.242
 0.025 0.654
 0.031 0.407
 0.032 0.243
 0.033 -0.166
 0.040 0.880
 0.056 0.019
 0.056 -0.121
 0.058 -0.265
 0.059 -0.586
 0.063 0.452
 0.065 0.676
 0.076 -0.582
 0.077 0.263
 0.091 0.693
 0.094 -0.233
 0.119 0.497
 0.123 0.385
 0.133 -1.089
 0.140 0.443
 0.149 0.427
 0.151 0.404
 0.162 1.012
 0.171 -0.138
 0.176 0.072
 0.186 0.743
 0.191 0.877
 0.199 0.458
 0.201 0.959
 0.202 -0.030
 0.204 0.493
 0.205 0.586
 0.210 0.583
 0.216 0.018
 0.216 0.545
 0.216 0.089
 0.219 -0.310
 0.231 0.397
 0.233 0.429
 0.237 0.085
 0.238 0.752
 0.240 1.151
 0.244 -0.216
 0.256 0.028
 0.269 0.246
 0.274 0.670
 0.275 0.255
 0.282 0.998
 0.288 -0.031
 0.290 0.644
 0.296 0.053
 0.297 0.733
 0.300 0.718
 0.302 0.244
 0.310 0.478
 0.321 0.880
 0.324 -0.148
 0.340 0.776
 0.349 0.675
 0.367 0.883
 0.373 0.398
 0.373 0.748
 0.377 0.842
 0.381 0.831
 0.381 0.692
 0.394 0.438
 0.395 0.680
 0.401 0.468
 0.401 1.038
 0.409 1.222
 0.419 1.393
 0.431 0.654
 0.440 0.694
 0.442 0.745
 0.446 0.542
 0.449 0.894
 0.454 0.688
 0.468 0.302
 0.471 0.655
 0.471 0.724
 0.484 0.968
 0.485 1.063
 0.486 0.900
 0.488 1.088
 0.493 1.032
 0.493 0.833
 0.494 0.998
 0.496 0.667
 0.498 0.742
 0.504 0.810
 0.505 0.941
 0.506 0.487
 0.512 0.606
 0.516 0.673
 0.519 0.884
 0.522 0.773
 0.523 0.686
 0.531 0.910
 0.540 0.915
 0.542 0.619
 0.549 0.851
 0.562 0.554
 0.567 0.511
 0.573 0.745
 0.581 0.632
 0.584 0.244
 0.586 0.708
 0.589 0.461
 0.601 0.793
 0.608 0.666
 0.610 0.399
 0.630 0.549
 0.637 0.298
 0.643 0.352
 0.647 0.662
 0.647 0.327
 0.648 0.519
 0.652 0.521
 0.657 0.074
 0.665 0.308
 0.667 0.237
 0.668 0.090
 0.668 0.368
 0.670 0.272
 0.671 0.467
 0.680 0.107
 0.683 0.304
 0.689 0.161
 0.689 0.244
 0.689 0.314
 0.690 0.379
 0.691 0.436
 0.698 0.212
 0.702 0.252
 0.703 0.132
 0.718 0.192
 0.718 0.189
 0.719 0.260
 0.731 0.126
 0.753 0.310
 0.768 0.313
 0.775 0.320
 0.779 0.345
 0.804 0.382
 0.805 0.347
 0.805 0.477
 0.817 0.370
 0.820 0.469
 0.822 0.525
 0.823 0.420
 0.826 0.366
 0.828 0.413
 0.830 0.407
 0.834 0.344
 0.839 0.339
 0.840 0.314
 0.845 0.300
 0.847 0.293
 0.848 0.265
 0.852 0.327
 0.857 0.218
 0.877 0.087
 0.880 0.066
 0.881 0.084
 0.885 0.011
 0.886 0.056
 0.886 0.020
 0.895 -0.048
 0.896 -0.049
 0.904 -0.101
 0.905 -0.071
 0.908 -0.071
 0.908 -0.118
 0.914 -0.092
 0.916 -0.113
 0.920 -0.116
 0.924 -0.120
 0.925 -0.098
 0.940 -0.088
 0.942 -0.072
 0.943 -0.092
 0.943 -0.078
 0.946 -0.075
 0.950 -0.056
 0.955 -0.018
 0.963 0.024
 0.970 0.066
 0.973 0.080
 0.975 0.089
 0.981 0.122
 0.982 0.128
 0.993 0.160
 0.996 0.164
 0.999 0.167
 ;

 title;

 proc sgplot data=Difficult;
    title "Scatter Plot of Y versus X";
    scatter x=x y=y ;
 run;

 title;

 proc iml;
    use difficult;
    /* read x and y from difficult into temp */
    read all var _num_ into temp;
    x = temp[,1];
    /* generate B-spline basis for a cubic spline
       with 4 evenly spaced internal knots in the x-range */
    bsp = bspline(x, 2, ., 4);
    Combined = temp || bsp;
    /* create a merged data set with x, y, and
       spline basis columns                    */
    create Combined var {x y col1 col2 col3 col4 col5 col6 col7};
    append from Combined;
 quit;

 proc ssm data=Combined opt(tech=dbldog);
    id x;
    /* parameters needed to define h(x) */
    parms v1-v7;
    /* defining h(x) */
    var = exp(v1*col1 + v2*col2 + v3*col3 + v4*col4
       + v5*col5 + v6*col6 + v7*col7);
    /* defining the polynomial spline trend */
    trend trend(ps(2));
    /* defining the observation noise with variance h(x) */
    irregular wn variance=var;
    model y = trend wn;
    output out=For pdv;
 run;

 proc sgplot data=For;
    title "Variable Bandwidth Smoothing Spline";
    band x=x lower=smoothed_lower_trend
       upper=smoothed_upper_trend ;
    series x=x y=smoothed_trend;
    scatter x=x y=y;
 run;

 proc sgplot data=For;
    title "Estimated Variance Function";
    series x=x y=var;
 run;