Try the commented and reproducible example below!
I use the `sf` package because is much faster than `rgdal` to open *Shapefiles*. Also, I used the great package `mapview` for interactive visualization in `R`. I transformed the points object to polygon's default projection before doing the intersection. 

**Note:** I downloaded the *Points.csv* file from Dropbox link first and saved it to working directory.


    # Load libraries ----------------------------------------------------------
    
    library('sf')
    library('mapview')
    
    # Download data -----------------------------------------------------------
    
    # Download polygon
    download.file("http://mapas.ine.pt/download/files/2011/nuts2/lisboa2011.zip", destfile = "lisboa2011.zip")
    unzip("lisboa2011.zip")
    
    # Load data ---------------------------------------------------------------
    
    polygons <- sf::read_sf("BGRI11_LISBOA.shp")
    points <- read.csv(file = "Points.csv", sep = ";")
    
    # Convert points df to sf object
    points <- st_as_sf(points, coords = c("lon", "lat"))
    st_crs(points) <- "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0" # assign CRS to points
    
    # Plot data
    mapview::mapview(polygons, zcol = "DTMN11") + mapview::mapview(points, color = "white", col.regions = "black")

[![plot][1]][1]    

    # Overlay Population with BGRI subsections --------------------------------
    
    # Transform points data
    pointsT <- st_transform(points, crs = st_crs(polygons))
    print(st_crs(pointsT))
    
    # Plot transformed data to check is ok
    mapview::mapview(polygons, zcol = "DTMN11") + mapview::mapview(pointsT, color = "white", col.regions = "black")
    
    # Intersection
    overlayOutput <- st_intersects(pointsT, polygons)
    
    # Build intersected polygons object
    intersectedPolys <- polygons[unlist(overlayOutput), ]
    
    mapview::mapview(polygons, zcol = "DTMN11") + mapview::mapview(pointsT, color = "white", col.regions = "black") + mapview::mapview(intersectedPolys, color = "black", col.regions = "black")

[![plot2][2]][2]

**Note:** the intersected polygons are colored in "black".

  [1]: https://i.sstatic.net/K5QQM.jpg
  [2]: https://i.sstatic.net/md0Nu.jpg