Hi, I am trying to generate RD plots with and without year FEs, but the graph doesn't look smoothed after the addition of FEs.

Code:
* Example generated by -dataex-. To install: ssc install dataex
clear
input int year float(running_var log_y uid)
2014  89.366  8.973331  1
2016  89.366  9.188673  1
2014  37.102  8.661576  2
2016  37.102  8.221179  2
2014  14.554  7.152594  3
2016  14.554  6.596597  3
2014   22.62  7.859138  4
2016   22.62  7.603361  4
2014  23.115  7.621384  5
2016  23.115  7.563481  5
2014  15.393  7.801198  6
2016  15.393  7.044147  6
2014   15.63  7.341237  7
2016   15.63  7.464042  7
2016  11.449  6.591065  8
2014  29.235  8.176641  9
2016  29.235   8.03603  9
2014  16.534  7.577722 10
2016  16.534   7.73528 10
2014  16.996  7.900858 11
2014 110.059  8.804719 12
2016  59.848  8.576594 13
2014   16.12  7.959042 14
2016   16.12  7.199755 14
2014 114.337  8.703652 15
2014  14.512  6.587711 16
2014  28.146  7.484158 17
2016  28.146  7.606704 17
2014   15.54  7.448243 18
2016   15.54  7.511616 18
2014 100.203  8.734136 19
2014  67.393  8.499454 20
2016  67.393  8.581483 20
2014  25.546  8.181071 21
2014  115.07  8.316578 22
2016  115.07    8.1924 22
2014  20.347  7.261263 23
2016  20.347  7.078819 23
2014  32.531   7.30103 24
2016  32.531  7.018701 24
2016  36.631  7.815246 25
2016  78.586  7.913814 26
2014  51.005 8.5136175 27
2016  51.005  8.789665 27
2014  46.133  7.859138 28
2016  46.133  7.309843 28
2014   19.01  9.022535 29
2016   19.01  8.146841 29
2014  73.146    8.1924 30
2016  73.146 8.6155815 30
2014   22.46  7.180985 31
2016   22.46  6.705008 31
2014  35.318  8.617797 32
2016  35.318  8.604086 32
2014  11.704  7.692583 33
2016  11.704  8.095413 33
2014  11.989  7.341632 34
2016  11.989  6.267172 34
2016  44.108  7.347915 35
2014  27.292  8.050651 36
2016  27.292  7.631241 36
2014 137.474  7.270213 37
2016 137.474  8.610766 37
2014  15.763  7.684396 38
2016  15.763  7.104487 38
2014  13.761  7.735838 39
2016  13.761  7.764027 39
2014  20.421  7.312177 40
2016  20.421   7.75151 40
2014 161.595  8.942698 41
2016 161.595  8.415491 41
2014  11.393  7.905742 42
2016  11.393  7.720407 42
2014  11.501  7.876102 43
2014  14.982   6.91169 44
2016  14.982  7.002598 44
2014  24.247  8.088526 45
2016  24.247  8.295765 45
2014  23.309  8.386142 46
2016  23.309 8.5717325 46
2014  10.428  9.830735 47
2016  10.428  8.745949 47
2014 129.653  8.197004 48
2016 129.653  8.308928 48
2014  14.022  6.466868 49
2016 108.738  7.598353 50
2014  72.212  7.382017 51
2016  17.383  6.932474 52
2014  25.646  7.209515 53
2016  25.646  6.599883 53
2014  18.232  7.209783 54
2016  18.232  7.653213 54
2016   11.03  7.223755 55
2014  28.753  8.414756 56
2016  28.753  8.177854 56
2014  33.928  8.362784 57
2016  33.928  8.072066 57
2014  39.279  8.005952 58
2016  39.279    8.2108 58
2014  90.547  7.536053 59
end
Code:
preserve
xi i.year
reg log_y c.running_var##c.running_var _Iyear* if running_var < 100, vce(cluster uid)
        predict yhat_left if running_var < 100, xb
        predict stdp_left if running_var < 100, stdp
       
        reg log_y c.running_var##c.running_var _Iyear* if running_var > 100 , vce(cluster uid)
        predict yhat_right if running_var > 100, xb
        predict stdp_right if running_var > 100, stdp

foreach side in left right {
        gen cipos_`side' = yhat_`side' + (1.96 * stdp_`side')
        gen cineg_`side' = yhat_`side' - (1.96 * stdp_`side')
    }

    keep running_var *_left *_right
    tempfile fit
    save `fit'

    * bin data
    restore, preserve
    
    qui: summ running_var
    egen bin = cut(running_var), at(0(10)200)
    replace bin = bin + 5
    collapse (mean) running_var log_y, by(bin)

    * --- Make plot
    append using `fit'
    twoway ///
        (line yhat_left running_var, sort(running_var) lcolor(black) lpat(solid) lwidth(medthick)) (line cipos_left running_var, sort(running_var) lcolor(black) lpat(shortdash) lwidth(thin)) (line cineg_left running_var, sort(running_var) lcolor(black) lpat(shortdash) lwidth(thin)) ///
        (line yhat_right running_var, sort(running_var) lcolor(black) lpat(solid) lwidth(medthick)) (line cipos_right running_var, sort(running_var) lcolor(black) lpat(shortdash) lwidth(thin)) (line cineg_right running_var, sort(running_var) lcolor(black) lpat(shortdash) lwidth(thin)) ///
        (scatter log_y running_var, xline(100) mcolor(black)), ///
        legend(off) graphregion(color(white)) xtitle("`xtitle'")  ytitle("`ytitle'") xline(100, lpattern(shortdash) lc(black)) ylab(, nogrid) note("`note'")
restore
Try running the regressions with and without "_Iyear*" term, and you will see the difference. What do I do to have the year fixed effects and have a smoothed graph as well?

Thanks!